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ABSTRACT 


This  is  the  Final  Report,  describing  all  the  tasks  accomplished  in  Phases 
A and  B of  Contract  No.  N00600-76-C-0790,  including  a listing  of  the 
computer  program  developed  under  the  present  contract  and  manual 
describing  input  and  output  data. 


In  order  to  incorporate  the  three-dimensional  effect  and  cascade  effect 
into  the  performance  prediction  of  the  super cavitating  propellers,  a two- 
dimensional  supercavitating  (2-D  s/c)  cascade  theory  and  a lifting  line 
theory  were  combined.  The  force  coefficients  obtained  from  the  2-D 
s/c  cascade  theory  will  account  for  the  cascade  effect  whereas  the  three- 
dimensionality  is  incorporated  in  terms  of  effective  flow  incidence  angles 
at  each  selected  spanwise  location  of  the  blade  for  the  2-D  program. 


An  inherent  difficulty  in  applying  the  2-D  s/c  cascade  theory  to  three- 
dimensional  flows  arises  due  to  the  existence  of  the  choking  condition 
but  was  overcome  by  correcting  the  effective  upstream  velocities  de- 
pending on  the  cavity  thickness.  Mathematical  formulation  combining 
the  2-D  s/c  cascade  theory  and  a lifting  line  theory  is  described.  The 
method  proposed  for  the  cavity  thickness  correction  is  explained,  fol- 
lowed by  numerical  procedures  to  solve  the  problem. 


Numerical  results  made  with  the  2-D  s/c  cascade  program  for  a s/c  cj;!  , • 
NSRDC  Model  3770  propeller  geometry  have  shown  a most  significant  

f 

i 

cavitating  cascade  effect.  These  results  seem  to  explain  very  well  the  • 


it 
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* discrepancy  existing  between  previous  experimental  data  and  design  data. 
The  propeller  characteristics  such  as  thrust,  torque  coefficients  and 
efficiency  have  been  calculated  and  compared  with  experimental  data, 

& having  provided  a good  correlation  over  a super cavitating  range  of  speed 

coefficient,  J. 

* However,  for  J' s beyond  the  above  range,  the  present  results  quickly 
deviate  from  the  experimental  data  because  a part  of  the  propeller  near 
the  hub  is  at  partially  cavitated  condition  to  which  the  present  theory  is 

* not  applicable. 
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INTRODUCTION 


The  development  of  a prediction  method  for  supercavitating  propeller 
performance  at  off-design  conditions  is  a difficult  task  due  to  an  additional 
complexity  of  the  cavity  flow  to  that  of  three-dimensional  propeller  con- 
figurations. 

Unlike  the  conventional  propellers  used  at  relatively  low  ship  speeds, 
supercavitating  propellers  are  expected  to  have  strong  cascade  effects. 
The  existence  of  blade  cavities  causes  blocking  or  choking  effects  on  the 
flow  passages  as  the  extent  of  cavity  becomes  large  both  in  length  and 
thickness.  Some  propeller  designers  already  pointed  out  the  importance 
of  the  cascade  effect  in  s/c  propeller  design  in  their  earlier  papers  such 
as  [l]  *. 

A similar  effect  was  also  found  important  even  for  subcavitating  propel- 
lers: the  paper  by  Kerwin  and  Leopold  [2]  showed  that  large  incidence 
angle  corrections  are  necessary  due  to  blade  thickness  effect  even  if  the 
thickness  is  small.  The  correction  becomes  particularly  significant  as 
the  thickness  ratio  to  the  blade  spacing  becomes  high,  i.  e. , near  the 
hub.  This  is  considered  to  be  exactly  the  same  blocking  effect  as  that 
for  the  cavity  flow.  It  is  now  evident  that  the  cascade  effect  of  blocking 
effect  must  be  correctly  incorporated  into  the  performance  prediction 
of  s/c  propellers. 

* 

Number  in  brackets  designates  Reference  at  end  of  paper 
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Although  linearized  s/c  cascade  theories  in  [3]  have  long  existed,  such 
theories  are  unable  to  accurately  predict  hydrodynamic  characteristics 
of  these  highly  nonlinear  s/c  cascade  flows.  It  was  not  until  just  re- 
cently that  a fully  nonlinear  2-D  s/c  cascade  theory  [4]  was  developed, 
greatly  facilitating  the  calculation  of  these  cascade  effects  and  providing 
a powerful  engineering  tool. 


The  method  of  solving  the  present  problem  is  a combination  of  the  2-D 

s/c  cascade  theory  with  a lifting  line  theory.  The  procedure  to  be  used 

is  stated  as  follows.  Specifying  all  physical  and  geometric  conditions  of 

the  s/c  propeller,  2-D  solutions  at  several  radial  or  spanwise  locations 

of  blades  will  be  obtained.  A difficult  question  however  arises  as  to  what 

effective  flow  incidence  angles,  ag , must  be  used  for  the  2-D  analysis. 

The  downwash  effect  in  propeller  flows  are  usually  so  strong  that  the 

geometric  flow  incidence  angles,  a (see  Figure  1 for  definition  and  also 

S 

Table  C2  for  actual  values  of  a for  the  3770  supercavitating  propeller), 

& 

are  completely  different  from  the  effective  incidence  angles,  ag. 


The  present  propeller  problem  is  very  much  similar  to  that  of  a single 
airfoil  of  finite  span  for  which  an  integral  equation  of  a lifting  line  theory 
must  be  solved.  The  result  determines  F,  the  distribution  of  circulation, 
or  equivalently  lift  over  the  blade  span  so  as  to  provide  a right  amount 
of  downwash  effect  everywhere  for  generating  the  above  circulation,  T. 
The  evaluation  of  the  downwash  angle,  a^*  in  this  case  is  most  simply 
made  by  a propeller  lifting  line  theory  but  in  a somewhat  complicated 
form.  The  effective  angle  of  flow  incidence,  ae»  is  then  obtained  by 
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II 


subtracting  a-  from  a . Applying  this  a to  the  lift  curves  calculated  by 
1 g e 

the  2-D  s/c  theory,  we  can  determine  the  circulation  distribution,  even- 
tually ending  up  with  an  integral  equation  for  ag  with  the  span  location  as 
a variable. 


A different  type  of  difficulty  arises  when  applying  a two-dimensional  flow 
approach  to  a three-dimensional  flow,  although  this  type  of  approach  has 
been  well  adopted  for  subcavitating  propeller  design.  Contrary  to 
the  supercavitating  propeller  problem,  the  same  method  for  subcavitating 
propellers  creates  no  serious  problems  in  determining  the  forces  at  any 
blade  location  for  any  given  effective  incidence  angle,  ae,  since  the  lift  and 
drag  forces  used  for  subcavitating  propellers  are  continuous  function  of 
a . However,  in  the  present  problem,  due  to  the  choking  condition  the 
force  curves  obtained  by  the  2-D  s/c  cascade  theory  are  discontinued 
right  at  that  point  (see  Figures  7(a)  to  (f)  for  choking  conditions  on  the 
lift  curves).  The  physical  meaning  of  this  is  explained  as  follows.  The 
cavity  length  and  thickness  increase  as  the  incidence  angle  increases, 
and  finally  the  cavity  extends  to  downstream  infinity  with  a maximum 
cavity  thickness.  This  blocks  or  chokes  the  flow  path  of  cascade.  It 
therefore  becomes  impossible  to  increase  the  total  mass  flow  going 
through  a cascade  beyond  that  point  at  the  choking  condition. 


This  type  of  2-D  choking  condition  never  occurs  in  the  three-dimensional 
(3-D)  flow  configuration  even  if  the  flow  cavitates  and  locally  chokes. 
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can  have  a similar  choking  condition  locally,  but  the  amount  of  flow  we  £ 

can  push  from  the  upstream  infinity  with  any  incidence  angle  is  unlimited 
since  any  mass  flow  in  excess  of  that  going  through  the  cascade  can  go 
around  the  corners  of  cascade  in  the  direction  of  span.  Thus,  the 
'effective'  flow  velocity  going  through  the  propeller  remains  almost 
constant  at  each  blade  radial  position.  The  terminal  values  depend  on 
the  cavity  thickness  but  not  depend  on  the  upstream  velocity.  As  long 
as  the  cascade  span  is  finite  this  phenomenon  holds  true.  However,  once 
the  span  extends  to  infinity,  going  back  to  a totally  two-dimensional  con- 
figuration, the  inherent  problem  mentioned  above  arises. 

As  a first  step  for  resolving  the  present  difficult  situation,  we  use  a 
simple,  intuitive  method  with  the  above  physical  picture  of  3-D  cavity 
flow  in  mind.  The  upstream  velocity  is  corrected  at  each  spanwise  lo- 
cation by  distributing  line  sources  in  cascade  configuration  whose  strengths 
are  determined  based  on  the  cavity  thickness.  The  effective  velocity 
obtained  with  this  method  is  always  smaller  than  that  of  the  original  flow 
so  that  the  cavitation  number  to  be  used  for  the  2-D  analysis  becomes 
larger  thus  being  able  to  avoid  the  choking  condition.  It  must  be  pointed 
out  that  the  correction  here  is  not  on  the  incidence  angle  as  downwash 
correction  but  on  the  upstream  velocity  or  equivalently  the  cavitation 
number. 

In  this  report,  we  present  a mathematical  formulation  which  combines 
the  2-D  s/c  cascade  theory  with  a lifting  line  theory  and  a method  for 
correcting  the  cavity  choking  effect,  followed  by  numerical  procedures 


4 


to  solve  the  problem.  Computed  2-D  s/c  cascade  results  for  a chosen 
NSRDC  Model  3770  propeller  geometry  are  then  presented,  showing  a 
remarkable  cascade  effect.  With  these  2-D  s/c  force  coefficients  used, 
the  propeller  performance  was  calculated  and  compared  with  the  ex- 
perimental data  [5]  . 

The  results  correlate  well  in  the  supercavitating  regions  but  quickly 
deviate  as  the  speed  coefficient,  J,  becomes  larger  due  to  the  appear- 
ance of  partial  cavitation  near  the  hub.  This  discrepancy  is  naturally 
expected  since  the  present  theory  is  only  applicable  to  the  fully  super- 
cavitating propellers. 

Originally,  the  present  work  had  been  planned  to  incorporate  the  results 
of  the  above  computations  into  a lifting  surface  theory  (see  Reference  [6] 
for  the  detailed  procedure  of  this  method).  However,  it  has  been  found 
that  the  present  method  combining  the  2-D  s/c  cascade  theory,  lifting 
lire  theory  and  cavity  thickness  correction  provides  accurate  results 
correlating  well  with  existing  experimental  data.  We  describe  several 
theoretical  backgrounds  why  the  method  accounts  for  all  supercavitating 
propeller  characteristics  as  follows; 

i)  By  having  used  the  results  of  2-D  s/c 

cascade  theory,  we  have  accounted  for 
a most  important  effect  of  super- 
cavitating propeller  flows,  i.  e.  the 
existence  of  the  cavity  in  cascade 
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geometry.  The  forces  on  the  s/c  cascade 
used  as  the  basis  for  propeller  calculations 
corrected  by  a lifting  line  theory  are  found 
to  be  much  smaller  than  those  of  a single 
foil  due  to  the  influence  of  the  low  pressure 
region  of  the  cavity  on  the  pressure  side 
of  an  adjacent  blade.  In  addition,  the 
blocking  effect  due  to  the  cavity  has  also 
been  incorporated  in  terms  of  a cavitation 
number  correction. 

ii)  The  finite  aspect  ratio  correction  to  s/c  pro- 
peller blades  is  usually  smaller  than  that  on 
fully  wetted  propeller  blades  (see  [8]  ) and 
also  limited  by  the  choking  condition.  There- 
fore the  lifting  line  theory  used  above  is  con- 
sidered fairly  accurate  as  this  has  been  proven 
by  a good  correlation  with  experimental  data. 

iii)  Although  a lifting  surface  theory  will  give  all 
boundary  value  corrections  including  both  the 
wetted  portion  of  the  blade  and  cavity  stream- 
lines as  mentioned  in  [6]  , no  information 
about  the  correction  for  upstream  flow  velocity 
is  obtained  by  the  present  method,  it  has  been 
proven  that  a correction  for  cavitation  number 
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is  a most  important  feature  in  applying  a 
2-D  s/c  cascade  theory  to  calculations  of 
the  propeller  performance. 

iv)  For  those  cases  in  which  the  cavitation  number  , 
a,  is  close  to  that  of  the  choking  condition,  the 
2-D  s/c  cascade  theory  itself  fails  to  converge 
in  the  numerical  iterative  procedure  as 
will  be  explained  later.  For  such  a's,  the  2-D 
theory  may  not  provide  a convergent  solution  for 
new  boundary  values  set  by  a lifting  surface  theory. 
In  the  present  approach,  however,  this  difficulty 
is  overcome  by  an  interpolation  scheme  as 
will  be  seen  later. 

Consequently,  we  believe  that  the  present  approach 
accurately  accounts  for  all  the  hydrodynamic 
effects  of  super cavitating  propellers  which  a lift- 
ing surface  theory  will  provide  and  furthermore 
that  the  former  is  superior  to  the  latter  from 
the  viewpoints  of  simplicity  in  concept  and 
economy  in  computation. 
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2. 


MATHEMATICAL  FORMULATIONS 


A two-dimensional  supercavitating  cascade  hydrodynamic  problem  has 
recently  been  solved  by  using  the  hodograph  variables  to  satisfy  the  exact 
boundary  conditions.  In  this  method  the  blade  and  cascade  geometry, 
the  upstream  flow  conditions  and  the  cavitation  number  are  specified. 

A system  of  five  nonlinear  functional  equations  involving  five  unknown 
solution  parameters  was  formulated  and  solved  numerically  using  a func- 
tional iterative  method  combined  with  Newton’s  method.  The  details  of 
the  theory  and  numerical  method  are  described  in  [4]  . 

In  order  to  incorporate  the  two-dimensional  (2-D)  cascade  theory  into 

the  analysis  of  supercavitating  (s/c)  propellers,  the  effective  angle  of 

incidence,  ag  (see  Figure  1),  must  be  determined.  The  geometric  flow 

incidence  a , which  is  determined  by  the  propeller  blade  pitch,  rotational 
§ 

speed  uu,  and  axial  flow  speed  V , is  typically  much  larger  than  a due  to 

the  strong  downwash  effects  generated  by  the  propeller  helical  vortex 

sheets.  For  example,  in  some  cases  of  s/c  propellers,  the  downwash 

angle  a-  = 3 - a can  be  as  high  as  ten  degrees  although  the  effective 

i g e 

angle  of  incidence  is  only  four  degrees.  It  has  become  clear  that  neglect- 
ing the  downwash  in  two-dimensional  cascade  calculations  can  result  in 
a solution  far  from  the  actual  propeller  flow  situation.  One  of  the  ideas 
in  capturing  the  three-dimensional  effect  is  to  incorporate  vortex  sin- 
gularities into  a lifting  line  theory. 


Preceding  page  blank 
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5{c 

With  the  geometry  of  cascade  and  propeller  blades,  U)  and  V specified,  % 

cL 


the  2-D  cascade  problem  can  be  solved  if  and  effective  upstream  velocity, 

V , are  assumed  known  at  each  radial  station  r since  a , 9 - and  a are 
6 6 1 

obtained: 

ai  = ag  ‘ ae 

(D 

0i  = 0 + ai 

(2) 

*\m*e 

a = r . 

kve 

(3) 

Five  equations  in  the  2-D  problem  are  now  rewritten  : 


fl  = 

-ae=0 

(Upstream  flow  angle  condition) 

(4) 

f2* 

Im  j.u(C^)  | + 2n  U^  = 0 

(Upstream  flow  velocity  condition) 

(5) 

f3a 

g3-v  o 

(Downstream  flow  angle  condition) 

(6) 

f , ■ 

s(l)-S  =0 

(Scaling  between  the  physical  and 

(7) 

4 

transform  planes) 

f5* 

gg-d  |sin^Xg  + ^-U2sin 

^2  + = 0 (Continuity  equation) 

(8) 

See  Nomenclature  and  also  the  Blade  Definition  Figure  1 for  the 
definition  of  each  symbols. 
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where  explicit  expressions  for  Re  uu( C ^ ) » Inn  uu(Cj)  > s( -1)  and  g$ 
are  given  in  Appendix  A (and  also  see  4 ).  It  is  noted  that  in  equations 
(4)  thru  (8)  ae's  simply  replace  in  equations  (7),  (8),  (9),  (15)  and  (16) 
in  4 . It  must  be  mentioned  that  the  upstream  velocity,  V , used  for 
these  2-D  calculations  i3  different  from  the  velocity  simply  composed  of 
the  axial  flow,  V , and  the  rotational  velocity,  aur,  as  is  shown  in  Figure 

cL 

1.  Due  to  the  three-dimensional  downwash  effect  and  the  cavity  blocking 
effect,  Vg  is  given  by  a following  equation: 


where  i 

Uj  = [((Dr  - wt)2  + (Va  + wa)2] 

wa»  w^  = propeller  induced  flow  velocities  in  the  axial  and 
tangential  directions 

= retarding  flow  velocity  due  to  the  cavity  blocking  effect. 

Before  describing  the  methods  of  determining  w , w.  and  V , first 
look  at  how  to  obtain  the  circulation,  T,  from  the  above  two-dimensional 
calculations,  which  will  be  used  in  a lifting-line  theory.  Taking  the 
control  volume  designated  by  ABCD  shown  in  Figure  2,  the  differences 
in  potential  between  the  points  A and  D,  and  B and  C are  calculated 
respectively  by 

A*DAa  vedsin(Y+  Ce' 

AcpCB  = U2  dsin  <Y  + a2^  * 
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thus  the  net  change  of  the  potential  in  this  control  volume,  that  is, 

F(x),  is  given  by 

r(x)  = &PDA  - AcpCB 

= Ve  dsin  (y  + ag)  - U2  dsin  (y  + a2) 

where  x = r/R. 

This  formula  holds  both  for  the  finite  and  infinite  cavity  cascade  flows, 
but  for  the  former  case  a simpler  form  is  obtained  by  using  a con- 
tinuity equation  between  the  upstream  and  downstream  flows,  i.  e. 

Vfi  d cos  (y  + aQ)  = U2  d cos  (y  + a2b 

sin(a  - cu) 

Hx}  = V d r-  , for  finite  cavity  flows. 

17  e cos  (y  + c^) 

F(x)  calculated  in  Equation  (10)  or  (11)  connects  the  2-D  results  with 
a three-dimensional  lifting-t:o.e  theory  to  find  the  propeller  induced 
velocities. 


The  induced  velocities  in  the  axial  and  tangential  directions  w and 
w^.for  the  case  where  the  blades  extend  from  the  hub  at  r = r^  to  the 
tip  r = R,  are  obtained  (see  [7]  for  detailed  derivations)  from: 


wa  1 

V"  = 2 

a 


dG(x') 

dx' 


x-x1 


ia(3.)dx- 


( 10  ) 


( ID 


( 12  ' 
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i and  i are  the  induction  factors  obtained  by  Lerb  [7]  and  detailed  ex- 

cL  C 

pressions  are  found  in  Appendix  B.  It  must  be  mentioned  that  for  s/c 
propellers  the  propeller  loading  is  expected  to  be  moderate  to  heavy, 
thus  the  downwash  effects  (12  ) and  (13  ) must  be  evaluated  by  taking 
into  account  the  deflection  of  the  vortex  sheet  location  behind  the  bound 
vorticies.  Lerb  [7]  showed  from  a discussion  of  energy  balance  be- 
tween the  propeller  disk  and  the  ultimate  wake  that  the  location  of  vor- 
tex sheets  should  be  on  a helical  surface  having  an  angle  3^  (instead  of 
9)  , which  is  a function  of  r.  In  the  present  calculations  of  w and  w 
we  use  0.  to  evaluate  i and  i . The  downwash  angle  a.  is  then  obtained 

1 cL  C * 

from  the  following  equation  to  an  accuracy  of  first  order  in  9.: 


. -1 

” an  ttx/J-w  /V 

t d 


a.  = 0. - 3 

i p i • 


( 16) 


13 


It  has  now  come  to  a point  of  how  we  incorporate  the  choking  or  cavity 
thickness  effect  into  the  problem,  which  relates  to  determining  V in 
Figure  1.  The  physical  picture  of  the  choking  phenomena  in  the  2-D  and 
3-D  flows  has  already  been  discussed  in  the  Introduction  of  this  report. 

A rigorous  treatment  of  this  type  of  problem  will  require  an  enormous 
effort  involving  complicated  mathematics,  although  it  must  be  done  some 
time  in  the  near  future.  Meanwhile,  we  use  a somewhat  more  intuitive 
method  as  a first  step  to  avoid  an  inherent  difficulty  in  applying  the 
results  of  2-D  s/c  cascade  flow  to  the  propeller  problem. 

In  order  to  represent  the  cavity  thickness,  a row  of  source  singularities 
of  strength  m are  placed  with  a distance  d in  a uniform  flow,  the  veloc- 
ity of  which  is  Uj  with  a stagger  angle,  V + a.Q,  as  depicted  in  Figure  3 . 

The  velocity  potential  of  the  flow  is  given  by 

W = Uj  Z e"l(Y+ae)  + m ln|sinh  )|,  (17) 

thus  the  velocity  potential  is  obtained; 

dW  TT  - 1( Y+  a ) , nan-  /.  . , „ . 

= Uj  e T e + -jj—  /tanh  (rrz/d).  ( 18  ) 

As  x ■*  +•,  the  x - component  of  the  velocity  changes  by  + mrr/d,  respec- 
tively. If  we  know  the  thickness  of  the  cavity,  d*  e,  the  strength  of  source, 
m,  is  calculated  by  using  the  continuity  equation 
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or 


{tJj  co„(v+ae>  - S|L}d=  {a2=0S(y+a2)  + T}'4-  de) 
Uj  cos(Y+ae)  - U2  COS  (Y+  a2)(  1 -e)  • 


niTT 

d 2 - e 


( 19) 


It  means  that  although  the  mass  flow,  Uj  cos  (y+  ae)*  d per  blade 
from  the  upstream  infinity, comes  into  the  cascade,  the  amount  of  mn 
is  rejected  to  go  through  the  blade  passage  due  to  the  existence  of  cavity. 
The  rejected  mass  flow,  miT,  should  go  normal  to  the  paper  plane,  in 
reality,  in  the  radial  direction  of  the  propeller.  Therefore  Vc  is  cal- 
culated from  (19)  by  taking  the  component  in  the  Uj  direction; 


V = 
c 


mrr 


d cos(y+ae) 


( 20  ) 


The  effective  upstream  flow  velocity  to  be  used  in  the  2-D  analysis  is 
now  obtained 


V = U.  - V„ 
e 1 c 


+ u 


cos  (y+  a2) 


2 cos(y+ 


a2  \ 
a ) / * 


1 - e 

2 - e 


( 21) 


where  e is  a function  of  a and  n,  obtained  from  the  results  of  the  2-D 

c 

computations.  Strictly  speaking,  the  present  method  is  only  valid  for 
the  infinite  cavity  flow  cases  in  which  the  cavity  is  fully  developed.  However 
even  for  the  finite  cavity  cases  it  is  considered  that  the  same  cavity 
blockage  evaluation  holds  true  by  taking  the  cavity  thickness  at  the  end 
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points  of  cavity  as  e in  Equations  (19)  and  (21). 

It  must  be  noted  here  that  the  correction  of  the  upstream  velocity  by  (21) 

changes  the  cavitation  number,  a,  for  which  e is  obtained  at  a . It  needs 

e 

an  iterative  scheme  to  satisfy  the  relationship  in  Equation  (21)  by  starting 
with  e for  a(U^)  and  ag  as  an  initial  step  and  then  finding  a new  e for  a 
new  a(V e)  where  V is  just  obtained  from  (21).  It  has  been  found  in  the 
actual  computations  that  the  convergence  of  the  iteration  is  rather  fast. 

The  problem  to  be  solved  is  now  fully  defined.  With  the  propeller  geome- 
try, V and  uj  specified,  one  can  determine  a circulation  distribution, 

a 

F(x),  in  such  a way  that  the  free  vortex  sheets  associated  with  the  T dis- 
tribution generate  a correct  amount  of  downwash  velocity  to  have  a 
sectional  blade  lift  equal  to  pU^T  where  Uffl  is  the  geometric  mean  velocity 
of  the  upstream  and  downstream  velocities  (see  Figure  4). 

It  is  immediately  seen  that  the  problem  is  completely  nonlinear  including 
integral  equations  and  thus  cannot  be  solved  explicitly.  Two  numerical 
iterative  methods  are  proposed  to  solve  this  type  of  situation  and  both 
procedures  used  here  will  be  explained  in  the  following  section. 
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3. 


NUMERICAL  PROCEDURES 


U 


Q 


& 


C; 


C 


Two  numerical  methods  for  solving  the  above  nonlinear  integral  equations 
are  proposed  and  have  been  tested  in  actual  computations  for  their  con- 
vergence. The  first  method  is  what  is  called  a substitutional  iterative 
method  and  the  second  one  is  Newton's  iterative  method  similar  to  that 
used  in  the  problem  of  three-dimensional  super cavitating  hydrofoils  [8]  . 

3.  1 Substitutional  Iterative  Method 

This  method  exactly  follows  the  steps  of  the  mathematical  formulation, 
the  flow  chart  being  shown  in  Figure  5. 

Assuming  the  effective  incidence  angles  a^te),  n=0,  at  each  spanwise 
location  one  can  find  downwash  angles  ct}  ;(x)  and  a cavitation  number 
<t(x)  from  equations  ( 1 ) and  ( 3 ).  The  solutions  of  the  two-dimensional 
s/c  cascade  problem  provide  the  deflected  flow  angle  at  down- 

stream infinity.  In  actual  computations  it  is  convenient  to  establish  a 
functional  relationship  of  as  a function  of  ae  and  a at  each  blade  sec- 
tion. Since  i-3  a smooth  function  of  and  g,  the  2-D  calculations 
for  several  values  of  ag  and  a* s will  be  sufficient  to  represent  by 
functionally  establishing  the  results  at  discrete  points.  By  doing  this 
one  can  save  a considerable  amount  of  computer  time  since  the  2-D 
computations  are  the  most  time  consuming  part  of  the  calculation.  If 
this  relation  is  not  established  at  the  beginning  of  the  computation  pro- 
cedure, the  2-D  program  must  be  run  for  each  iterative  loop.  This  can 
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be  seen  in  the  flow  chart.  Figure  5,  where  the  returned  loop  v/ill  now 


go  back  to  the  2-D  calculation  box  instead  of  the  a^-box.  In  some  cases 
in  which  the  stagger  angle  and  blade  solidity  are  large,  the  2-D  s/c 
cascade  program  becomes  numerically  unstable  as  was  reported  in  [4]  . 
This  problem,  however,  is  overcome  if  the  2-D  features  are  completely 
calculated  at  the  initial  stage  of  the  numerical  procedure. 


The  induced  velocity,  V c,  due  to  the  existence  of  cavity  and  thus  the 
effective  flow  velocity,  V , are  obtained  by  a small  iterative  procedure 
in  Equation  (21).  The  sectional  circulation  distribution  T(x)  is  then  ob- 
tained by  Equation  (10)  or  (11),  thereby  enabling  us  to  calculate  w , w. 

3>  t 

and  a^n+1^(x).  The  values  of  a^°^(x)  first  assumed  are  now  checked  to 

determine  that  they  are  corrected.  If  not,  with  a new  a^n+^(x)  and 

a^n+l),  we  proceed  to  the  next  iteration  until  a convergent  solution  is 

( n) 

obtained.  In  each  iteration,  0.  , starting  with  an  assumed  value 

/ \ 1 

= 0 + a^°  y,must  be  calculated  and  a new  value  of  0^n^  must  be  used 

in  calculations  of  w and  w. . It  must  also  be  noted  that  the  cavitation 

a t 

number  a based  on  V is  used  for  the  first  iteration  but  a based  on  V 


is  used  from  the  second  iteration  on. 


e 


e 


If  the  test  for  the  convergence  of  solution  parameters,  for  example,  a^, 
is  pa  ,sed,  we  proceed  to  calculate  the  propeller  characteristics  such 
as  thrust,  power  coefficients  and  efficiency. 

When  the  method  was  applied  to  the  present  problem,  we  found  that  con- 
verged solutions  were  obtained  only  if  assumed  starting  values  of  a ^ 

6 


is 


were  close  to  the  actual  solutions.  It  is  for  this  reason  that  a second 
method  using  Newton's  technique  is  proposed  for  seeking  a better 
convergence. 

3.  2 Newton's  Iterative  Method 

We  incorporate  Newton's  method  into  the  nonlinear  integral  equations 
for  improving  the  convergence  of  iteration.  This  requires  a new  arrange- 
ment of  the  problem  in  order  to  identify  the  solution  parameters. 


From  Equations  (16),  (B-16)  and  (B-17), 


f s tan  (8  (x)  - a (x) 
— \ g e 


k 

-{l  + TTi-^  = Gmhm*  (»(x))}  - 0 


(22) 


and  from  Equations  (14)  and  (B-ll), 

Ji,  n*>  »,<*).  oa) 

S *L,  Gmsinm®(*)  - = 0 

m=  1 a 


( 23  ) 


Choosing  discrete  control  points  in  the  radial  direction  of  the  blade  for 
which  the  computations  will  be  made,  say  x = 0.4  to  0.9  by  0.1  increment, 
we  have  six  independent  equations  in  (22)  so  that  the  same  number  of 
Gj^'s  are  chosen  for  the  solution  parameters,  in  this  case  k=  6.  Since 
all  other  quantities  in  equations  (22)  and  (23)  are  known  except  for 
ae(x)'s,  they  are  naturally  chosen  as  another  six  solution  parameters. 
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called  aek«  We  now  have  2k  solution  parameters  for  a system  of 
nonlineai  integral  equations  having  an  order  of  2k. 

Rewriting  these  equations  and  parameters  symbolically  by 

E = & s) 

(Gk'  aek)  » 

we  can  describe  the  above  set  of  equations  as  follows: 

F(x)  = 0 

thus  Newton's  iterative  loop  is  established  by 

J (X(n))  . (x(n+ - X(n))  = - F(x(n))  (24) 

where  J is  a Jacobian  matrix  whose  component  is  given  by 

J-  aFj/aXj.  (25) 

In  the  present  case  each  component  of  J is  either  analytically  or  numer- 

— 

ically  calculated; 


A(x.)  6ik/cOS  ■ aek}  ; 1 = 1 ~ 6' 

j = 7 ~ 12,  k = j - 6 


ST=\  = sin{jcp(xk)}:  1-  7~  12,  j»  1-6,  k=i-6 


35"^  " 2nR’\r-  * ^ r (xk'  aek'  ae);  1 = 7 ~ 12« 


j s 7 ~ 12,  ksi.6,  Z-j-6 


J 


where  all  partial  derivatives  are  analytically  calculated  except  for  T 
for  which  a finite  difference  method  is  used. 

Iterative  numerical  procedures  for  this  case  shown  in  Figure  6 are  very 
similar  to  those  of  the  first  method  shown  in  Figure  5.  Our  experience 
in  using  this  Newton's  method  for  the  present  problem  indicated  rather 
slow  but  steady  convergences  for  almost  all  cases.  It  has  also  been 
found  that  the  method  is  much  less  sensitive  to  the  initial  starting  values 
of  solution  parameters. 
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4. 


CALCULATIONS  OF  THRUST,  TORQUE  COEFFICIENTS  AND 


EFFICIENCY 


In  the  cascade  flow  the  lift  force  acting  on  the  blades  is  known  to  be 
normal  to  a geometric  mean  angle  a.  (see  References  [4]  and  [9])  which 
is  depicted  in  Figure  4: 


where 


«.  = <=<’»'1j?g;(veC05a<» + u2cosa2)} 

O.  = * [ V«2  + U22  + 2ve  U2  009  (ae  ‘ “2) 


V is  taken  to  be  unity  in  the  2-D  calculations  and  from  Figure  1 


V /V  is  calculated  from 
e a 


U,  Vrt  + Vr  / wt\  ! 
ir  = -v^  = (cot3-va)cos0i  • 

Thus,  a sectional  thrust  is  obtained: 


dT  = g|pUwT  cos  ^i+ag-a.^-  (DCB+Df)3in^i+ae"a»)|  dr 
where  D and  Dt  are  pressure  drag  on  the  cavitating  blade  parallel  to 

CO  I 

the  direction  of  U^  and  friction  drag  on  the  propeller  blade; 

D«  = CD»  * C 2 P U0D 


Df  = Cf . c \ p Ua 
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#r  cos  I 0.  + ae  - CLm  j - (D^  + Df)  sin  (0. 


= g 2ttR2  p Va2 


1 

Jj  y-  (x')  G<x')  cos  ^ +’ae  - 


(X1)  /U 


^ *R-  \-J  (CI>+Cf)sin(0 


i + ae  ’ °J  } 


The  thrust  coefficient  C^,  is  obtained  by-  normalization: 
t r i UJX’>  / 

/-•  - J-  I i ® \ / o i — 


CT  ■ 7 t — r = 4sf  {“V“ 

ipV  2 TTR2  J ( Va 

a *h 


G(x')  cos 


(3i  + «e  - ac) 


i c(x,)  l^Y  ,c 

4n  R \ Va  / (C] 

sin  (0i  + ae  - a„)  | dx'  . 


<CD.+  Cf>  X 


Similarly  the  power  coefficient  is  calculated  as: 


dP  = ru)  dF  dr 


= rtug  | pU^r  sin  ^0.  + afi  - aB  j + (D^  + Df)  cos  (0j  + a0  - a*  ) | dr 


P = g 2ttu)R3  o V 2 

3> 


f (U0D(X')  / \ 

J x'  ) “V G(x,)  sin  (8i  + ae  ‘ a- J 

xi_  * a \ / 


L * Va 


, c(x')  / Ujx')\2  > 

h?  ) (CD-+Cf)  cos(0i  + ae-aco)|dx' 


C = 


P ipV/  TTR4 


= Jx' {“V"- G(x')  sin  (9i+ae-^)) 
xh 

1 c(x,)/U-  \2 

+ TFT  r \ va  / (CD-  + Cf)  x 
cosO.  + ag-a.jjdx' 


COS (3; 


r 


where 


It  must  be  emphasized  that  U is  calculated  as  a nondimensional  number 
in  the  2-D  cascade  theory  problem,  referring  it  to  V . However,  in 
Equations  (33)  to  (36),  need  to  be  absolute  values.  They  must  be  multiplied 
by  Ve  in  these  calculations.  Similarly,  T calculated  in  the  2-D  problem 
by  equation  (10)  must  use  ’ d*  which  has  a dimension  since  again  in  2-D 
calculations  d is  normalized  by  the  chord  length  c.  The  propeller  ef- 


g ficiency  r|  is  finally  calculated  as: 


CT 

’-Cp  • 

( 37) 

• 

Another  definition  of  thrust  and  torque  coefficients,  using  symbols  K^, 
and  ^Q,is  given  by 

t 

Kt  = T/pn2D4 

( 38  ) 

Kq  = Q/pn2D5  , 

( 39  ) 

thus  the  relations  between  these  numbers  and  C_,  and  C are  obtained 

G T P 


TC  — C* 

KT  " ~T  CT 

(40) 

KQ  = T5-  Cp 

(41  ) 

and 

kt  j 

(42) 

C 
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An  alternative  way  to  obtain  these  coefficients  is  that  instead  of  using  T 
and  we  can  directly  use  and  which  are  normal  and  parallel 
to  the  direction  of  the  upstream  velocity,  V . C-,  and  C are  now  ex- 

® 1 P 

pressed  by  the  following  formulae: 

} /V  (x')  \2  ( ) 

CT  = 2 J s°l(3e' ) I — — J x'  <CL  cos  0.-CD  sin8. -Cfsin0g>dx' 

Xh 

Cp  = -|  J'sol(x')  ^ ^ x'2  |cL  sin0.+  CD  cos  0.  + Cf  cos  @gjdx’ 

*h 

where  sol(x')  is  a solidity  of  the  blade  at  x. 


( 43  ) 


( 44  ) 
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NUMERICAL  RESULTS  OF  THE  TWO-DIMENSIONAL  SUPER- 
CAVITATING  CASCADE 


As  shown  in  the  flow  charts  of  Figures  4 and  5 the  first  step  in  the  numerical 
procedure  of  the  present  method  is  the  calculation  of  the  two-dimensional 
hydrodynamic  characteristics  of  super cavitating  cascade. 

» 

It  is  intended  that  results  of  trial  computation  will  be  compared  with 
existing  experimental  data  [5]  and  those  being  currently  obtained  at 
David  Taylor  Naval  Ship  Research  and  Development  Center  (DWTNSRDC). 

The  geometry  of  a super  cavitating  propeller  Model  3770  designed  on 
the  basis  of  the  method  developed  by  DWTNSRDC  [ij  has  therefore  been 
chosen.  The  profiles  of  the  blades  were  designed  based  on  a Tulin- 
Burkart  two-term  camber  with  an  additional  camber  to  account  for  a 
lifting  surface  correction  [lo]  . Appendix  C describes  the  equations  of 

^ the  two-term  camber  with  some  representative  coordinates  and  K,  cor- 

rection factors,  including  other  hydrodynamic  design  and  geometric  pa- 
rameters. 


In  order  to  cover  a complete  matrix  of  possible  effective  incidence 
angles  ag  and  local  cavitation  numbers  oq  with  J in  the  2-D  computations, 
the  following  procedure  is  used. 


C 


The  design  cavitation  number  of  Model  3770,  based  on  the  ship  speed 
V ,is  chosen  to  be  0.617.  First  of  all,  the  local  cavitation  number 

f 2 2)  2 

based  on  V(=<(cur)  + V > ) can  be  calculated  at  each  radial  location: 


Pi -Pc 
v ' lpV2 


= a 


Va 


|l  + (nx/J)2}  • 


(45 


av's  calculated  this  way  are  listed  in  Table  C2  of  Appendix  C. 


Strictly  speaking,  however,  the  cavitation  number  to  be  used  for  2-D 
calculations  must  be  based  on  the  effective  upstream  flow  velocity  Vg 
so  that 


where  V /U,  and  V /V  are  calculated  from  Equations  (32)  and  (21), 
alee 

respectively.  As  a matter  of  fact,  these  0"e ’ s have  been  used  in  the  pres- 
ent propeller  computations. 


From  the  flow  angles  3,  blade  setting  angles  3 and  speed  coefficient 

J,  the  geometric  incidences  angles  a , can  be  calculated  to  estimate 

the  initial  values  for  a . Since  0 is  shown  in  the  Table  Cl  of  Appendix 

e g 

C and  0 is  calculated  from: 


0 - tan"  *(  J/rrx), 


(46) 


a is  easily  obtained  and  is  tabulated  in  Table  C3. 
g 
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It  is  seen  from  Table  C3  that  the  maximum  and  minimum  values  of 

are  21. 11  degrees  and  2.  06  degrees  at  x=  0.  3,  J=  0.  3 and  x=  0.  9,  J=  0.  7, 

respectively.  The  range  of  cavitation  number  based  on  V in  Table  C2 

is  found  to  be  0.  0069  to  0.  2194.  Based  on  these  values  it  was  decided 

to  calculate  the  2-D  s/c  cascade  characteristics  at  four  different  incidence 

angles,  2,  3,  4 and  6 degrees  with  a cavitation  number  ranging  from  the 

choking  condition  to  about  0.08.  For  any  other  combination  of  an  incidence 

angle  and  a which  will  arise  in  the  iterative  procedure,  the  2-D  flow 

characteristics  will  be  extrapolated  or  interpolated  analytically.  It  is 

noted  that  the  maximum  value  of  a.Q,  i.  e.  6 degrees  does  not  seem 

to  cover  a value  of  a of  21.  11  degrees.  However,  the  downwash  effect 

§ 

near  the  hub  is  so  large  that  the  effective  angle  of  attack  will  be  near  or 
within  6 degrees.  It  is  also  obvious  that  no  supercavitation  occurs  at 
a = .2194. 


Figures  7(a)  thru  7(f)  show  the  calculated  lifts  normal  to  the  upstream 

flow  as  functions  of  cavitation  number  a at  normalized  radius  locations, 

e 

x=0.4,  0.5,  0.6,  0.7,  0.8  and  0.9.  The  2-D  calculations  were  left  out 
for  the  point  at  x=  0.  3 because  the  cavitation  number  is  too  large  and  the 
solidity  is  too  high  to  obtain  convergent  solutions  in  the  2-D  computations. 
In  addition,  the  propeller  performance  can  be  accurately  calculated  without 
the  information  at  that  point  by  an  interpolation  scheme  as  will 
be  seen  later. 


Two  different  computer  programs  (see  [4]  ) were  used,  one  for  the  chok- 
ing condition  at  which  the  cavity  extends  downstream  to  infinity  and  the 


r. 
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other  for  the  finite  cavity  case.  In  these  figures  we  see  that  a signifi- 
cant cascade  effect  occurs  in  cavity  flows.  In  Figure  7(a),  for  example, 
where  the  solidity  is  small,  0.244,  near  the  tip  (at  x=  0.9)  with  a stag- 
ger angle  of  74  degrees,  it  is  seen  that  the  lift  coefficient  increases 
as  the  incidence  angle  increases.  This  phenomenon  is  quite  similar  to 
that  observed  in  single  lifting  foil  cases.  It  means  that  the  solidity  of 
0.  244  at  this  location  with  y=  74°  is  yet  too  small  to  see  much  of  a cas- 
cade effect  and  thus  the  flow  is  similar  to  that  of  a single  foil  except  that 
the  choking  phenomenon  appears.  However,  at  x=  0.  8 where  the  solidity 
becomes  slightly  larger,  0.365,  with  a stagger  angle  of  72.4  degrees, 

the  lift  coefficient  CT  at  a = 6°  loses  its  value  as  ct  becomes  small  (see 

Li  e 

Figure  7(b)),  until  finally  its  value  becomes  even  smaller  than  that  ob- 
tained at  a0  = 4°,  3°  and  2°.  The  reason  why  this  occurs  at  smaller  a' s 
is  obvious:  the  smaller  the  cavitation  number,  the  longer  and  thicker 
is  the  cavity  (see  Figures  10(a)  thru  10(f)),  so  that  the  cavity  boundary 
with  a low  cavity  pressure  is  close  to  the  pressure  side  of  the  neighbor- 
ing blade,  causing  a loss  in  lift.  It  is  also  seen  that  the  cavitation  number 
CT  at  which  this  change  in  CT  occurs  in  Figures  7(a)-(f)  checks  quite 
well  with  the  value  of  ag  at  which  the  length  of  cavity  starts  extending 
to  infinity  (choking  conditions)  as  shown  in  Figures  10.  One  can  also 
observe  a similar  behavior  in  for  ae  = 4°,  occurring  here  at  a smaller 
CT  than  for  the  ae=  6°  case.  This  trend  becomes  even  stronger  (see 
Figures  7(c)  thru  7(f))  since  the  solidity  becomes  larger  increasing  from 
0.479  to  0.912.  In  Figure  7(f)  where  x=  0.4  and  the  largest  solidity 
occurs,  the  relation  between  the  lift  and  incidence  angle  completely 
flips  over  for  a range  of  CTe's,  i.  e.  the  lift  is  the  highest  at  the  lowest 
incidence  angle. 
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This  peculiar  behavior  for  cavitating  cascade  flow  observed  above  never 
happens  in  the  cases  of  single  lifting  foils.  Physically  it  can  be  under- 
stood and  explained  as  follows.  When  the  solidity  becomes  large  and 
thus  the  blades  are  more  closely  packed  in  the  cascade  configuration, 
the  adjacent  blades  are  strongly  affected  by  the  existence  of  cavity  thus 
causing  significant  hydrodynamic  effects.  This  effect  becomes  stronger 
as  the  cavity  becomes  thicker  and  longer  or  as  the  flow  incidence  angles 
become  larger  and  the  cavitation  number  becomes  smaller  as  has  been 
seen  above. 

To  our  knowledge  the  above  highly  nonlinear  cascade  effect  have  never 
been  incorporated  into  super  cavitating  propeller  design.  If  the  lift 
curves  of  single  super  cavitating  foils  are  used  for  such  designs,  large 
discrepancies  between  the  design  and  experimental  data  are  to  be  ex- 
pected. For  example,  increasing  flow  incidence  angles  or, equivalently, 
increasing  blade  camber  at  a radial  location  of  the  blade  having  a re- 
latively large  solidity  essentially  decreases  the  sectional  lift.  This 
results  in  a smaller  thrust  coefficient  in  experiments  than  expected 
by  design.  A totally  opposite  situation  must  sometimes  be  taken;  blade 
angles  and  camber  must  be  decreased  to  increase  lift  depending  on  the 
solidity  and  cavitation  number.  This  may  be  one  of  the  reasons  why  the 
experimental  thrust  a^id  torque  coefficients  of  NSRDC  Model  3770  were 
found  far  short  of  the  design  values  based  on  single  foil  predictions. 
More  detailed  discussions  about  this  point  will  be  made  in  the 
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next  section. 
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It  is  interesting  to  compare  the  lift  coefficients  obtained  in  the  present 
nonlinear  cascade  theory  and  those  of  super cavitating  single  foils.  The 
latter  values  at  ag  = 0 are  easily  computed  from  the  design  lift  coef- 
ficients, correction  factors  listed  in  Table  Cl  and  angles  of  attack: 


• K + TTOC./ & 


(47) 


where  the  subscripts  S and  0 in  (47)  designate  'single  foil'  and  'zero 
cavitation  number',  respectively.  CLS0  calculated  based  on  Equation 
(47)  are  plotted  in  each  Figure  7(a)  through  7(f).  It  is  seen  that  a well 
known  approximation  for  finite  cavity  length  by  a correction  factor  (1  + a), 
commonly  used  for  a single  foil  flow, cannot  be  applied  to  the  s/c  flow  in 
the  cascade  configuration  whatsoever.  It  is  also  noted  that  C^^q's  are 
much  larger  than  those  values  extrapolated  from  the  linear  portions  of 
curves,  again  showing  a remarkable  super  cavitating  cascade  effect 
on  lifting  forces. 

It  is  also  seer,  that  the  choking  conditions  marked  in  Figures  7(a)  through 
7(f)  vary  to  a great  degree  depending  on  the  solidity.  With  small  solidity 
and  a small  incidence  angle,  the  choking  flow  does  not  occur  until  ag  be- 
comes fairly  small,  say  0.007  (see  Figure  7(a)),  while  a of  0.041 
is  enough  to  cause  the  same  condition  for  a large  solidity  and  a large 
angle  of  attack  (see  Figure  7(f)).  This  behavior  is  also  attributable  to 
the  increasing  cascade  blockage  effects  with  increase  in  solidity  and  in- 
cidence angle. 
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Figures  8(a)  through  (£)  show  drag  coefficients  parallel  to  the  upstream 
flow  direction,  each  corresponding  to  Figure  7(a)  through  (f),  respectively. 
It  is  seen  that  these  drag  forces  also  exhibit  a trend  similar  to  that  of 
the  lift  coefficient. 


The  lift  and  drag  coefficients  shown  in  these  figures  will  be  used  later 
in  propeller  analysis  to  calculate  thrust  and  torque  coefficients  by  using 
the  formulae  in  Equations  (43)  and  (44).  The  information  which  now  con- 
nect the  2-D  theory  to  the  propeller  lifting  line  theory  are  those  about 
the  circulation,  P.  Using  the  computer  outputs  of  the  2-D  s/c  cascade 
theory,  in  particular  the  flow  deflection  angle,  a>  and  U^/V ' , we  calcu- 
lated T's  and  plotted  them  in  Figures  9(a)  to  (f).  Equations  (10)  and  (11) 
are  used  for  P's  of  the  infinite  and  finite  cavity  cases,  respectively. 

The  numbers  read  out  from  these  figures  are  used  as  input  data  to  a 
computer  program  for  s/c  propellers. 

Finally,  Figures  11(a)  through  11(f)  show  lift-to-drag  ratios  (L/D)  as  a 
function  of  lift.  It  is  interesting  to  see  that  the  values  of  L/D  are  less 
sensitive  to  cavitation  number  (or  in  the  figures)  and  solidity  as  long 
as  incidence  angles  and  other  parameters  remain  constant.  This  indicates 
that  a blade  section  having  a good  L/D  value  as  a single  foil  with  an  infinitely 
long  cavity  is  also  guaranteed  to  have  a good  L/D  at  finite  cavity  lengths 
in  a cascade  configuration.  This  fact  also  seems  to  explain  the  good 
correlation,  obtained  in  the  efficiency  of  the  3770  propeller,  between 
experimental  data  and  design  data,  while  the  thrust  and  torque  coef- 
ficients are  way  off  as  already  discussed. 
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The  numerical  computations  presently  carried  out  with  the  two-dimensional 
super cavitating  cascade  programs  are  a most  time-consuming  and  difficult 
part  during  the  present  analysis.  In  particular,  such  flow  configurations 
as  having  large  solidity  with  large  stagger  angles  and  large  incidence  an- 
gles cause  numerical  instabilities  in  the  functional  iterative  procedure 
as  was  pointed  out  in  [4]  . In  the  present  case,  for  example  at  x=  0.4 

and  a = 6°,  where  the  solidity  is  0.  912,  the  next  value  of  a to  the  choking 
e 

case  which  could  be  calculated  was  o = 0.  068.  For  any  cavitation  num- 

e 

ber  between  these  two  points,  the  numerical  procedure  failed  to  ob- 
tain a convergent  solution.  However,  other  points,  where  a > 

0.  068,  were  calculated  without  any  difficulties  and  these  points  have  been 
smoothly  connected  by  a curve  as  shown  in  Figure  7(f).  The  error  in- 
curred in  this  way  is  not  considered  to  be  significant  in  the  present 
analysis.  For  all  other  cases,  convergent  solutions  were  obtained  at 
almost  equally  spaced  values  of  Execution  time  to  obtain  a conver- 

gent solution  at  each  data  point  was  about  150  seconds  with  CDC  6600 
and  about  40  seconds  with  CDC  7600.  About  sixty  data  points  were  com- 
puted to  generate  the  present  2-D  s/c  cascade  data,  so  that  a total  of 
2400  seconds  of  computer  time  was  used  with  the  CDC  7600  or  equiva- 
lently 9000  seconds  with  the  CDC  6600. 


34 


6 


CALCULATIONS  OF  PROPELLER  PERFORMANCE  FOR 
NSRDC  MODEL  ;>770  SUPERCAVITATING  PROPELLER 


Based  on  the  mathematical  formulation  and  numerical  procedures  de- 
scribed in  the  precedent  sections,  a computer  program  has  been  written 
for  calculating  K,p,  Kq  and  (thrust,  torque  coefficients  and  efficiency) 
of  supercavitating  propellers.  (A  complete  listing  of  the  computer  pro- 
gram and  input-output  data  manual  are  given  in  Appendix  D). 

The  2-D  s/c  cascade  data  for  the  3770  propeller  have  been  already  pre- 
pared for  propeller  analysis.  By  taking  five  discrete  radial  points  on 
the  blade,  i.  e.  , x=0.4,  0.6,  0.7,  0.8  and  0.9,  the  propeller  hydro- 
dynamic  characteristics  have  been  calculated.  The  cavitation  number, 
CTya  and  speed  coefficient,  J,  of  the  3770  propeller  at  the  design  point 
are  0.617  and  0.44,  respectively  (see  Table  1 of  Reference  5)  . First 
of  all  we  calculated  K,p,  Kq  and  ri  at  this  point  and  show  the  results  in 
Table  2 in  comparison  with  design  and  experimental  data  taken  from 
Reference  5 (also  see  Figure  12).  It  is  clearly  seen  that  the  present 
method  predicts  them  well,  in  particular,  thrust  coefficient  K^,,  and 
efficiency  r\  particularly,  i.  e.  within  4 percent  of  the  experimental  data. 
There  exists  a large  discrepancy  in  K^,  and  Kq  between  the  design  data 
and  experiment  which  use  about  15  and  11  percent,  respectively,  although 
the  r|  there  is  close  to  others.  The  reason  for  this  discrepancy  has  al- 
ready been  explained  in  the  previous  sections;  the  data  basis  for  the  pres- 
ent method  depends  on  the  supercavitating  cascade  theory  whereas  the 
design  method  depended  on  an  infinite  cavity,  single  foil  theory. 

I 
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This  point  will  become  even  more  clear  when  we  compare  some  of  blade 
sectional  characteristics  between  the  two  methods.  Four  different  pa- 
rameters are  shown  in  Table  3 for  comparison.  These  include  local 
cavitation  number,  effective  flow  incidence  angle,  downwash  angle  and 
lift  coefficients.  The  local  cavitation  numbers  of  the  present  method 
shown  in  the  table  are  those  corrected  on  the  basis  of  the  cavity  thick- 
ness data  from  the  2-D  cascade  calculations  (see  Equation  (21)),  whereas 

2 2 - 

those  of  all  design  methods  are  simply  based  on  V = ((uur)  + V )2. 

a 

The  discrepancies  shown  in  all  these  parameters  of  Table  3 are  quite 
large  and  it  is  again  considered  that  they  are  attributable  to  the  differ- 
ences in  the  force  characteristics  used  by  the  two  different  methods  as 
was  mentioned  before.  Among  others,  it  is  seen  that  one  of  the  most 
significant  differences  exists  in  downwash  a.  angles  and  thus  effective 
flow  incidence  angles  aQ  , the  latter  in  the  design  method  are  all  about 
2 degrees  whereas  the  former  range  from  6.8  degrees  to  4.6  degrees. 

The  local  lift  coefficients  neai  the  hub  predicted  by  the  present  methods 
are  larger  than  those  of  the  design  method  but  become  smaller  as  one 
proceeds  toward  the  propeller  tip.  A question  now  arises  why  the  thrust 
and  power  coefficients  predicted  here  are  smaller  than  those  of  the  design 
method.  It  is  simply  answered  that  lift  forces  near  the  tip  are  more  con- 
tributory to  K^,  and  Kq  due  to  the  smaller  pitch  and  larger  radius  as 
will  be  readily  seen  from  Equations  (34)  and  (36). 

It  is  interesting  to  investigate  the  boundary  between  the  supercavitating 
and  partial  cavitating  regimes  of  propeller  operation.  According  to  the 
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2-D  s/c  cascade  data  regarding  the  length  of  cavity,  shown  in  Figures 

10(a)  to  (f),  it  reads  that  the  propeller  blade  at  x=  0.4  is  partially  cavi- 

tated  or  fully  wetted  if  ae=  2.  06°  and  <7=  . 0656  obtained  by  the  design 

method  are  true  (see  Figure  10(f)).  With  a = 6.  8°  and  0=  . 0798  at 

e 

x=  0.4  predicted  by  the  present  method,  the  length  of  cavity  is  about 
1.  1 chord  length  which  means  that  the  propeller  barely  stays  in  a super- 
cavitating  regime.  Unfortunately,  no  clear  photographic  evidence  of 
3770  is  available  for  us  to  check  this  feature.  In  any  case,  it  is  strongly 
recommended  that  not  only  force  data  but  also  photographic  data  at  each 
measurement  point  be  taken  for  any  future  supercavitating  propeller  ex- 
periments. The  former  will  justify  the  overall  accuracy  of  a prediction 
method  used  whereas  the  latter  will  play  an  important  role  in  verifying 
a local  flow  phenomena. 


Figure  12  shows  K,p,  Kq  and  p over  a range  of  speed  coefficient,  J,  for 
the  design  cavitation  number  (7ya  = .617.  As  J is  taken  clcse  to  about 
0.  5,  the  calculated  a and  a at  x=  0.4  become  4.  5°  and  . 097.  Again 
from  Figure  10(f),  we  can  see  that  the  blade  near  the  hub  (at  x=  0.4) 
under  the  above  conditions  is  in  a partially  cavitating  region  where  the 
present  supercavitating  propeller  theory  becomes  invalid.  It  seems 
this  reason  that  all  curves  in  Figure  12  start  showing  deflection  around 
such  a J. 


In  Figure  13  we  compare  the  present  results  with  experimental  data 
for  <7ya  = . 500.  It  is  seen  that  the  correlation  for  K,p,  Kq  and  p at 
such  J of  0.44  to  . 5 is  excellent  but  that  the  discrepancy  starts  growing 
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as  J becomes  outside  the  above  range  of  J.  The  same  reason  as  above 
explains  the  discrepancy  for  large  J's  at  which  a part  of  the  propeller 
operates  at  partial  cavitating  condition. 

The  reason  for  the  discrepancy  for  smaller  J (say,  less  than  0.44)  is 
not  known  at  this  point.  Physically,  a decrease  of  J increases  the  ef- 
fective flow  angle  a and  decreases  local  cavitation  number  a , leading 
more  and  more  to  a local  choking  condition  of  supercavitating  propeller. 

From  certain  values  of  a and  a on,  neither  of  these  values  cannot  change. 

e e 

It  means  that  the  lift  and  drag  coefficients  should  reach  constant  values 

thus  C,p  and  C in  Equations  (43)  and  (44)  also  approach  constant  values. 

2 3 

However,  K^,  and  Kq  are  proportional  to  J and  J respectively  (see 
Equations  (40)  and  (41)  so  that  these  values  theoretically  decreases  as  J 
decreases.  Experimental  data  in  Figure  13  do  not  show  any  of  this  type 
of  trend  at  least  before  J=  0.4  whereas  the  results  of  the  present  theory 
clearly  demonstrate  the  above  theoretical  argument. 

tt  is  also  considered  that  the  discrepancy  in  K,p,  Kq  and  r|  for  smaller 
J may  be  suggestive  even  for  incidence  angle  correction  due  to  the  choking 
condition  of  the  2-D  cascade  theory  in  addition  to  the  a-correction.  For 
determining  which  argument  is  true,  further  comparisons  of  the  present 
theory  with  existing  and  possibly  new  experimental  data  from  super- 
cavitating propellers  are  urgent. 

In  order  to  accurately  predict  propeller  performance  for  a range  of  larger 
J,  say  0.  51  in  this  case,  the  2-D  cascade  data  based  on  partial  cavitating 
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conditions  must  be  used.  Although  a linearized  theory  for  partial  cavi- 
tating  cascade  is  available,  it  cannot  be  applied  for  the  propeller  analysis 
in  which  the  flow  there  is  highly  nonlinear.  Development  of  a similar 
nonlinear  theory  to  that  for  supercavitating  cascade  4 is  now  in  order. 

Figure  13  also  shows  results  before  we  applied  the  correction  method 

to  cavitation  numbers  (by  dotted  lines).  When  J was  set  at  a slightly 

smaller  value  than  0.46,  in  this  case,  a combination  of  such  a and  a 

e 

put  us  beyond  the  choking  condition  where  there  exists  no  2-D  cascade 
data.  By  forcing  as  to  find  T by  extrapolation  (which  seems  irrelavant, 
first  of  all)  we  managed  to  determine  K^,  Kq  and  p,  which  are  plotted 
in  Figure  13.  These  values  are  of  course  not  valid.  Without  using  either 
cavitation  number  correction  as  does  the  present  method  or  possibly 
incidence  angle  correction,  it  is  impossible  to  avoid  this  inherent  2-D 
choking  problem  in  three-dimensional  flow  applications.  It  is  clearly 
seen  that  an  adequate  treatment  of  this  problem  is  most  crucial  in  the 
entire  supercavitating  propeller  analysis  as  long  as  a cascade  theory  is 
intended  to  be  used. 

For  the  above  computations,  we  have  used  Newton's  iterative  procedures 
which  showed  a slow  but  steady  convergence.  For  e«.  ch  run,  it  took  about 
30  seconds  for  a CDC  7600. 
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7.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  present  method  has  incorporated  a 2-D  s/c  cascade  theory  into  a 
propeller  lifting  line  theory  with  downwash  angle  correction  and  cavita- 
tion number  correction  made. 

The  results  obtained  from  the  theory  have  shown  an  excellent  correla- 
tion with  experimental  data  within  a certain  range  of  speed  coefficient,  J, 
whereas  some  discrepancy  exists  outside  that  range.  It  is  quite  clear 
that  the  discrepancy  for  larger,  J,  is  attributable  to  the  appearance  of 
the  partial  cavitating  flow  near  the  hub  so  that  the  present  theory  naturally 
becomes  invalid.  That  for  smaller,  J,  cannot  be  clearly  explained  at 
this  stage  as  mentioned  before. 

Furthermore,  it  has  been  found  that  there  exists  a significantly  large 
difference  in  local  flow  conditions  between  the  present  theory  based  on 
cascade  data  and  design  theory  based  on  single  foil  data.  These  include 
downwash  flow  angles,  effective  flow  angles,  cavitation  numbers  and 
lift  coefficients. 

In  order  to  clarify  some  of  the  uncertainties  having  just  arisen,  the  fol- 
lowing specific  recommendations  for  future  research  works  to  be  car- 
ried out  are  made. 

1)  More  comparisons  between  the  theory  and 

experiments  for  different  types  of  propeller 


Preceding  page  blank 
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o 1 

configurations  such  as  those  in  11  are  necessary 
for  verifying  the  accuracy  of  the  present  prediction  theory. 

2)  Slightly  beyond  a range  of  the  largest  ef- 
ficiency  or  optimum  design  point,  there 
exists  a partially  cavitating  flow  regime. 

We  always  have  a possibility  that  s/c 
propellers  sometimes  operate  at  such  a 
regime.  In  order  to  cover  a full  range 
of  performance  prediction  for  s/c  pro- 
pellers, it  is  advisable  to  develop  a non- 
linear partially  cavitating  cascade  theory 

3 I 

and  to  incorporate  the  data  from  such  a 
theory  into  the  present  method. 

3)  Although  the  two-dimensional  supercavitating 
cascade  theory  [4]  used  here  was  compared 
well  with  experimental  data  before,  such 
a comparison  was  very  much  limited  due 
to  the  lack  of  experimental  data.  It  is 
recommeded  that  more  experiments  be 
conducted,  in  particular  for  a range  of 
high  solidity  and  high  stagger  angle. 

| 

4)  A more  rigorous  evaluation  of  the  effects 
of  propeller  cavity  thickness  on  cavitation 
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number  and  possibly  on  local  flow  angle 
may  be  necessary. 


Finally,  it  is  concluded  that  cascade  and  three-dimensional  effect  of 
supercavitating  propellers  plays  a most  crucial  role  in  their  hydrodynamic 
performance.  In  cooperation  of  this  effect  into  both  prediction  and  design 
methods  for  such  propellers  is  unavoidable. 
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Figure  3 Representation  of  Cavity  Thickness  by 
Distribution  of  Source  Singularities  in 
a Cascade  Row. 
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Figure  4 
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Geometric  Mean  Angle  q^.  Force  Normal 
to  Geometric  Mean  Velocity  U®  in  Cascade 
Flow  and  Thrust  and  Torque  Components. 
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(start) 


2-P  S/C  Cascade  Program 

• Propeller  Geometry#  »<«7l  specified 

• a^n*  assumed  (n»0)  -*  0jB^»  9^  1 **  2) 

• found  (a*  0)  (Eq.  3) 


Obtain  (Geometry,  a^B*. 


whore  c in  a desired  value 
vtor  convergence  .. 


Calculate  K^.  K p fc  ^ 
(Eqs.37.  38  it  39) 


STOP 


Figure  5 Flow  Chart  of  Substitution  Procedures  Iterative 
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STOP 


Figure  6 Flow  Chart  of  Newton's  Iterative  Procedures 
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Figure  7(a)  Lift  Coefficient  C-^  Normal  to  the  Upstream 

Velocity  vs.  Cavitation  Number  a for  Incidence 
Angles  ae=  2°,  3°,  4°  and  6°  at  the  Blade  Span- 
wise  Location  of  x=  0.  9 where  the  Solidity  is 
0.244  and  geometric  stagger  angle  y is  74.03°. 

A O □ are  CL  Values  of  ag  = 2°,  4°,  6° 
Respectively  Calculated  from  a Linearized  Theory 
for  a Single  Foil  (Equation  47). 


51 


C Figure  7(e)  The  Same  as  Figure  7(a)  Except  That 

x=  0.  5 where  the  Solidity  is  0.  728 
and  y is  63.  94°. 
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Figure  9(c)  F/dVgVS,  ae Corresponding  to  Figure  7(c) 
(x=  0.7). 


Figure  10(b)  The  Same  as  Figure  10(a)  Except  That 
x = 0.  8. 


70 


J 


CAVITY  LENGTH  /CHORD 


0 .02  .04  .06  .08  .10 


Figure  10(e)  The  Same  as  Figure  10(a)  Except  That 
x=  0.  5. 
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Figure  11(a)  Lift-to-Drag  Ratio  L/D  vs.  for  Incidence 
Angles  ae  = 2°,  4°  and  6°  at  x = 0.  9. 
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THRUST  COEFFICIENT  TORQUE  COEFFICIENT  EFFICENCY 


C J = .44 

Design  poinf 


SPEED  COEFFICIENT  J 


Figure  12  Performance  Prediction  for  3770  at  Design  <7ya=  .617. 
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TABLE  2 

Comparison  of  K^,,  Kq  and  ri  between  the  Present  Result, 

Design  Data  and  Experimental  Data  for  3770 

Supercavitating  Propeller  at  Design  Point, 

a.r  = . 617  and  J = .440. 
v a 


a 

Design  Data 
(Ref.  5 ) 

Experimental  Data 
(Ref.  5 ) 

Present 

Results 

kt 

. 1004 

.085 

.0819 

kq 

.0130 

.0115 

.0106 

T i%) 

54.  1 

52.0 

54.0 
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APPENDIX  A 


Two  Dimensional  Nonlinear  Super cavitating 
Cascade  Theory  (Reference  [4]  ) 


The  detailed  mathematical  formulation  is  described  in  [4],  However 
some  important  features  needed  for  the  present  calculations  are  re- 
peated here  for  convenience. 

Rejx(Cj)|  , Imjx(£i)|  , g^,  s(-l)  and  g5  in  equations  (4)  through  (8)  are 
given  by: 


uKC1)  = + i ImjtttCj)} 


= V(C 
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-h 


28(5') 
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2TTiJ  v (?'+i)(5'-b) 
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where 


Cj  = A exp  ji(rr/2-6)|  , 
g3  = ^ in(\  l+a/U2)*  j.2n: 
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h^l,  a,b,c,U2(a2)J  =[exp|" 
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k(?,A)  = S 


§ cos  6 


(g-A  sin  6)^+  (A  cos  5)^ 
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and 


5 = a + Y » 
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(A-6) 


d I -ifi  Crc  16  £i~cl 
§5  = 2??ie  in  C^+e 


(A-7) 


In  these  equations  a,b,c  are  § -coordinates  in  the  mapped  plane  (see 
Figure  A-2),  A is  a parameter  associated  with  the  cascade  mapping  func- 
tion and  U,  is  the  velocity  at  downstream  infinity  which  is  related  to  a2 

u 

through  a continuity  equation: 


A-2 


1 


¥ 


cos(CL  + Y) 

T T — ^ 

u 2 cos(a2+  Y ) 


where  is  taken  to  be  unity  in  the  2-D  calculations  and  all  of  these  values 
are  solution  parameters  to  be  determined  by  Equations  (4)  through  (8) 
shown  in  the  text.  The  potential  plane  from,  which,  the  C -plane  is  mapped 
and  the  definition  of  s and  3 are  shown  in  Figures  A1  and  A3. 


APPENDIX  B 


Calculations  of  Induced  Velocities  w_^  and  by 
a Lifting  Line  Theory 
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The  following  calculation  method  of  the  induced  velocities  w and  w is 

di  c 

based  on  the  work  of  Lerbs  [7].  The  equations  for  w and  w are  written 

cL  C 

again: 


w (x) 
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1 dG(x^_^  0.,x',x)dx' 

2 dx'  x-x'  a i 


(B-l) 


wfc(x) 


where 


i 

if  dG(x')  _1_  i ,g 
2 x-x-  V V 

•V 


, x)  dx' 


a 


( B-2  ) 


o 


G(x)=  2TTRV 


w = 


x 

S x'tan  3 


x-  tan  3 


X < X- 


x > x' 


( B-3  ) 


(B-4) 


W* 


*(f  -1)B2- 

-*(jr  -1)<l+Bl>> 


x<  x' 


X > X- 


(B-5) 


B-l 


MW f* 


B 


'Azfif  ln  (1+  /i.zj. 


(B-6) 


1,2=  ‘(^ 


1+y2.  Vl+y' 


/ 2 (Vi+y1  +D(\i+y  -l) 


( B-7 ) 


^ tan  t 


( B-8  ) 


^ x'tan 


(B-9) 


where  Nicholson's  asymptotic  formulae  have  been  applied  in  obtaining 

B and  B-  from  the  original  integrals  of  vortex  sheets. 

1 2 


By  introducing  cp  for  a change  of  variables: 


x = Z (1  + xh^“  2 (1-Xh^ 


COSC0  , 


( B-10  ) 


thus  x = 


*h 


and  1 correspond  to  sp  = 0 and  it. 


We  also  write  G(x)  and  i in  Fourier  sine  and  cosine  series,  respectively: 
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It  must  be  noted  that,  at  cp=0  and  cp=TT: 
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where  L'  Hospital's  rule  has  been  used. 
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APPENDIX  C 

Geometric  and  Hydrodynamic  Configurations 
Super cavitating  Propeller  Model  TMB(NSRDC) 


of 

3770 


C.-i 


Sectional  profiles  of  supercavitating  propeller  Model  TMB(NSRDC)  3770 
are  made  up  by  Tulin-Barkart  two-terms  camber  sections  at  zero  cavi- 
tation number,  modified  with  a lifting  surface  correction  factor.  The 
equations  of  such  profiles  are  given  by: 


£=  ^-k(4(S)+|  d)3/2-  4<i!2 

c 5rr  (3'c  3 ' c c 


where 


K 


design  lift  coefficient 
chord  length 

correction  factor  of  a lifting  line  surface  theory  • 


In  Table  Cl  the  y-coordinates  of  the  blade  pressure  sides  at  several 
radial  locations  are  shown  with  appropriate  values  of  and  K.  Also 
shown  in  the  same  Table  are  the  pitch-to -diameter  ratio,  P/D,  blade 
setting  angle,  3^,  geometric  stagger  angle  y and  solidity,  including  the 
number  of  propeller  blades,  diameter  of  propeller  and  hub-to-tip  dia- 
meter ratio. 


Table  C2  shows  the  geometric  flow  angles  S (see  Figure  1)  and  the  geo- 


metric flow  incidence  angles  ctg  at  various  radius  locations  for  various 

speed  coefficients,  J,  It  is  noted  that  the  sum  of  B and  or.  is  equal  to 

§ 

the  geometric  blade  setting  angles  8 at  any  location  for  all  J’s. 

8 


The  local  cavitation  numbers  based  on  V(  = j(a'r)^+ V^2)  are  calculated  by 
equation  (45  ),  for  various  J's  and  are  shown  in  Table  C3.  This  table 
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roughly  shows  a range  of  the  cavitation  number  over  which  the  super- 
cavitating  propeller  hydrodynamics  are  to  be  calculated. 
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Table  C-l  Geometric  and  Hydrodynamic  Configurations 

of  Super cavitating  Propeller  Model  TMB(NSRDC)  3770 
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APPENDIX  D 

Computer  Program  Listing  and 
Input  and  Output  Data  Setup 
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INTRODUCTION 


The  computer  program  called  'SCSCREW'  (listed  below)  calculates  the 
hydrodynamic  characteristics  of  super cavitating  propellers  with  sectional 
two-dimensional  s/c  cascade  data  given  as  input  data.  Therefore,  the 
computer  program  developed  in  [4]  must  be  used  to  generate  these  2-D 
s/c  cascade  data  prior  to  the  use  of  'SCSCREW'.  The  method  of  preparing 
these  input  data  will  be  explained  later. 

In  what  follows  we  describe  the  structure  of  the  program  'SCSCREW' 
including  functions  of  subroutines,  input  data  set-up  and  type  of  output 
data  obtained  as  the  result  of  calculations. 
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2.  STRUCTURE  OF  SCSCREW 


SCSCREW  consists  of  a main  program  and  several  subroutines,  brief 
descriptions  of  which  will  be  given  as  follows: 


1)  MAIN  PROGRAM  SCSCREW 

• Specify  the  dimensions  for  data. 

• Read  input  data. 

• Exercise  Newton's  iterative  procedure  (see  Figure  6). 

• Calculate  Cj,,  C^,  Kj,,  Kq  and  T|  with  and  without  drag  forces. 

• Calculate  local  flow  conditions  including  downwash  flow  angle, 
effective  incidence  angle,  cavitation  number, lift  and  drag 
coefficients. 


2)  FUNCTION  CLCD  (I,  S,  B,  ILD) 

• Interpolate  lift,  drag  and  circulation  with  input  data  passed 
on  to  this  program  through  common  statement. 


I:  Index  for  radial  or  spanwise  position  on  blade 

S:  Cavitation  number  for  which  CLCD  to  be  calculated 

B:  Flow  incidence  angle  for  which  CLCD  to  be  calculated 

ICLCD:  Control  Index 


ICLCD  = 0 for  lift 


ICLCD  = 1 for  drag 
ICLCD  = 2 for  circulation 


Preceding  page  blank 
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3)  SUBROUTINE  FITCAV  (MM1,  SE,  EE,  Ml,  II,  ACC) 

• Curve-fitting  for  cavity  thickness  by  polynomials. 

MM1:  Index  for  radial  or  spanwise  position  on  blade 

SE:  Cavitation  number  at  which  cavity  thickness  data  are 

available 

EE:  Cavity  thickness  data 

Ml:  Number  of  flow  incidence  angles  for  which  cavity 

thickness  data  are  available 

II:  Number  of  cavitation  numbers  for  which  cavity 

thickness  data  are  available 

ACC:  Coefficients  of  polynomials  fitted  for  cavitation  number 

as  a function  of  cavity  thickness,  calculated  in  FITCAV. 

4)  SUBROUTINE  CAVINO  (NG,  SC,  AE,  ACC,  A2,  II,  SCN,  EE,  Ml) 

• Correct  a local  cavitation  number  for  the  effect  of  cavity 
thickness  (see  Equation's  (17)  through  (21)). 

NG: 

SC: 

AE: 

ACC: 

A2: 

II: 

SCN: 


Index  for  radial  or  spanwise  position  on  blade 
Cavitation  number  before  correction 
Flow  incidence  angle 

Coefficients  for  polynomials  with  which  cavitation 
number  is  fitted  as  a function  of  cavity  thickness 
in  Subroutine  FITCAV 

Flow  incidence  angles  at  which  cavity  thickness  data 
are  available 

Number  of  cavitation  numbers  for  which  cavity 
thickness  data  are  available 

Corrected  cavitation  number 


■ 

i 

} 
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EE: 


Cavity  thickness  data 

Number  of  flow  incidence  angles  for  which  cavity 
thickness  data  are  available 


r 

re. 

[C 

f 

;C 

c 


t> 


Ml: 

5)  SUBROUTINE  FANC  (F,  FE,  SC,  ACC,  NG,  KI,  II) 

• Provides  a functional  relation  for  cavitation  number  as  a 
function  of  cavity  thickness; 

f h (SC  • U12/Ve2)  - (Polynomials  in  FITCAV) 

where  Uj/Ve  is  given  in  Equation  (21)  of  the  text. 

F:  f given  above 

EE:  Cavity  thickness  (e  in  Equation  (21)) 

SC:  Cavitation  number 

ACC:  Coefficients  of  polynomial  in  FITCAV 

NG:  Index  for  a radial  position  on  the  blade 

KI:  Index  for  an  incidence  angle 

II:  Dummy  index  (to  be  neglected) 

6)  SUBROUTINE  MOSEC  (A,  B,  ER1,  ER2,  X,  J,  NG,  KI,  SC,  ACC,  II) 

• Find  a root  for  f(x)  = 0 where  x must  lie  between  A and  B and 
f( A)  > 0,  f(B)  < 0. 

A,B:  A root  of  f(x)  = 0 lies  between  A and  B 

ER1,  ER2:  Accuracy  control  valuables  with  which  |xrea^  - x |<ER1 

and  |f(xreal)  - f(x)  [<  ER2 

x:  A root  for  f(x)  = 0 found  in  this  subroutine 

J:  Number  of  iterations  executed  in  MOSEC 

NG,KI,  SC,  ACC,  II:  The  same  as  those  in  FANC. 
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7)  SUBROUTINE  DETERM  (A,  N,  D) 


— 


• Calculate  determinant  of  a matrix  A of  rank  N. 

A:  Matrix  input,  requiring  dimension  (50,  2N+3) 

N:  Rank  of  the  matrix 

D:  Calculated  determinant  of  A 


8)  SUBROUTINE  LSQUAR  (DATA,  NUMBER,  N,  A,  CHISQ,  XM) 

• Least  square  fitting  for  DATA  with  polynomials  of  order  N, 

DATA:  DATA(  1 , NUMBER)  = x 
DATA(2,  NUMBER)  ■ y 
DATA(3,  NUMBER)  = Error  in  data 

NUMBER:  Number  of  input  data 

N:  Order  of  polynomials 

A:  Coefficients  of  polynomials 

CHISQ:  Chi-square  error  to  be  specified 

XM:  Dimension  of  (20,  43)  needed,  but  neglect  data  in 

XM  (not  used  here) 

9)  FUNCTION  FALF(X),  FGAM(X)  & FBET(X) 

• Calculate  a,  8,  y at  x in  Filon' s integration  formula  (see 
Equations  (25.4,47)  to  (25.4.  57)  in  "Handbook  of  Mathematical 
Functions",  National  Bureau  of  Standard). 

10)  SUBROUTINE  SPLINE  (X,  Y,  DY,  S2,  S3,  T,  SS,  SSI,  SS2,  L1M,  N,  Cl). 

• Cubic  spline  curve  fitting  for  Y(X)  and  evaluate  SS(T) 

Y(X):  Dimensional  Data  of  order  N 

T:  Points  for  which  Y to  be  evaluated 

SS:  Evaluated  data  a T 

OTHER  PARAMETERS:  Disregard  (not  used) 


j 


i 
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3.  INPUT  DATA  SET-UP 


The  following  describes  set-ups  for  input  data  cards.  Typical  input  data 
L.  are  also  listed  at  the  end  of  the  computer  program  listing. 


Card  No.  Symbol  Description  FORMAT 

1 NGAUS  Number  of  points  for  8110 

Gauss  Quadrature 

2 T(I)  Nondimens ional  positive  4F20.  10 

half  coordinates  for  Gauss 
Quadrature 

3 W(I)  Weighting  factors  for  4F20.  10 

Gauss  Quadrature 

4 EP  Increment  for  a finite  8F10.  5 

difference  method  to  ob- 
tain partial  differentials 
in  Newton's  method 

5 MAXIT  Maximum  iteration  number  8110 

in  Newton's  method 

MM  Number  of  discrete  con- 

trol points  on  the  propeller 
blade 

MF,  NF  Number  of  terms  used  for 

Filon  integration  method 

NFILON  Number  of  increments  in 

Filon  integration  method 

IOLD  If  not  equal  to  0,  VIVA  of 

old  calculations  are  fed  in 
as  input  data.  If  0,  it  is 
approximately  calculated 
in  the  program  (VIVA  s U^/Va) 

IWRITE  Number  of  the  last  iterations 

for  which  output  printing  is 
made. 
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Card  No. 

Symbol 

Description 

FORMAT 

6 

DIA 

VASHIP 

Propeller  diameter  in  inches. 

Propeller  axial  speed  in 
feet/sec  (=  V ) 

a 

8F10.  5 

7 

XJJ 

zz 

SIGVA 

XH 

XXM 

Advance  coefficient  (J) 
Number  of  blades  (g) 
Cavitation  number  based  on 
Va  <=  °Va» 

Propeller  hub  diameter/ 
propeller  tip  diameter  (=  x^) 

Weighting  factor  in  iterative 
procedure 

8F10.  5 

8 

XX(I) 

Nondimens ional  radial  or 
spanwise  position  (=x)  where 
I from  2 to  MM-1 

8F10.  5 

9 

BETAG(I) 

Geometric  blade  angle  (see 
Figure  1)  in  degree  ( = 8) 

5 

8F10.  5 

10 

ALFE(I) 

Effective  flow  incidence  angle 

in  degree  ( = a ) 
e 

8F10.  5 

11 

SOLI(I) 

Solidity  (=  SOL) 

8F10.  5 

12 

VIVA(I)  U1  / Va  data  °nly  if  IOLD  * 0 

— Repeat  for  EXX  from  2 to  MM-1 

8F10.  5 

13 

XXX 

Nondimens  ional  radial  position  F10.3 
(=x)  for  lift,  drag  and  circulation 
data 

MANGLE 

Number  of  incidence  angles 
(=CL  ) for  which  the  data  are 
available 

110 

ISIG 

Number  of  points  for  a for 
which  the  data  are  available 

no 

SIGMIN(IXX) 

Not  used,  disregard 

F10.  5 
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Repeat  MANG1  times 


¥ 


14 

ANGl(IANG) 

Incidence  angle  at  which 
the  data  are  read-in. 

8F10.  5 

15 

SIGD(IXX,  IANG, 
I) 

Oe  for  which  the  data  are 
read-in  at  IXX  & IANG 

8F10.  5 

16 

CL1(IXX,  IANG, 

Data  for  lift  coefficients 

8F10.  5 

I) 

it 

n 

17 

CDL  (IXX,  IANG, 
I) 

Data  for  drag  coefficients 
(=  CD> 

8F10.  5 

18 

GGGI( IXX,  IANG, 

Data  for  circulation 

8F10.  5 

(=r/dv  ) 

e 


Repeat  for  IX X = 2 to  MM-1 

19  XXX 

MANG1  ) The  same  as  before 
ISIGl 


8F10.  5 


20 


21 


22 


-Repeat  for  IANG  = 1 to  MANGl 
ANG2(IANG) 


Incidence  angle  at  which  8F10.  5 

cavity  thickness  data  are 

read-in 


SIGE(  IXX,  IANG,  Cavitation  number  at  which  8F10.  5 
I)  cavity  thickness  data  are 

read -in 

EE(IXX,  IANG,  I)  Cavity  thickness  at  IXX  and  8F10.  5 
IANG 
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TYPICAL  OUTPUT  DATA 


Typical  output  data  are  also  listed  at  the  end  of  the  program  listing. 
Most  of  them  are  self-explanatory.  Those  not  explained  in  output  data 
are  described  as  follows: 


C 


€ 


G 


C 


XN(I)  : solution  parameters 

XN(I)  | = ae  , 1=1,  MM 

| = G , I = MM  + 1,  2 x MM 
1 m 

F(I):  Residues  of  each  function  in  Equation  (22) 

P(I,  J):  Partial  derivatives  of  Jacobian  J in  Equation  (25) 

ALFG(I):  a 

o 

BET  AI(I):  0. 

ALFKI):  a. 

VIVA(I):  Ul/Va 

SIGV(I):  local  cavitation  number 


5 COMPUTER  PROGRAM  LISTING  OF  'SCSCREW' 


PROGRAM  Sc3CReW( INPUT, OUTPUT,  TAPE 5* INPUT * TAPE6«0UTPUT ) 
t LlFTTNG  LINE  PROPfLLCR  THEURT  in  COMBINATION  pith  2-0  3/C  CASCAOE  THEORY. 
L PpOGpAMMLD  BT  0.  FljRUYA,  5-27-76. 

DIMENSION  hETA(IO),  HETAIdO),  ALFE (10)*  VIVAtlO),  GGAM(IO) 
nlHENSION  GGG(99),SnH(10),ALFINF(10),XX(10).ALFEO(10) 

DIMENSION  PH1(I2),  T(100),  W(IOO),  SUUT ( 1 50 ) , GM ( 1 09 ) 

DIMENSION  3U5o)*<JFTD(100),OFPUUOO)  ,UFCP(100),ALFI(10) 

DIMENSION  cLS(l0),CO3(10)*ALFG(10),FTO(m,FPOU0) 

9 1 PENSION  y ALF (10)*  XUET(IO),  SINAI  10  • 150) * S INT ( 1 0 , 1 50 ) , BET AG( I 0 ) 
DIMENSION  hMA(iG*100)>  HMT(10*100),  WAVA(IO),  WTVA(IO),  SIGV(tO) 
DIMENSION  ijtVI(lo),  FCT ( 12) , FCPI12),  XCTI 12) * OFCT ( 1 00 ) *UI V A(10 ) 
DIMENSION  nr(l50)*S2(150)*S3(I50),SSI(150),SS2(150) 

DIMENSION  XIAC8,150),XIT(B,I50),ALFINO(10) 

dimension  GMUl20),Pf50,20),0(SO,20),F(20) ,XN(20) ,BGA(20) , TBGI20) 
DIMENSION  DrTE(20),PXJ(20),SlGVU(20),VIVAO(20)*CDLIM(10,10) 

DIMENSION  ANG2(l0),STGE(10,5,5).EE(t0,5*5),ACC(10»5»5) 

DIMENSION  sHt(tO) 

COMMON  PAI,  COnV,  CrJNV  I , SIGMJNI 1 0 ) 

COMMON  SIGOUO,5,15),  CL1 ( 1 0 . 5, 15) » ANG1 (5) »C0t (10,5* 15) 

COMMON  MANGLE,  IS IG* GGGI (10*5,15) 

pA 1*3. 141 592654 

CONV*PAI/tflO. 

CONVI*180./PAl 


* 


C 


L OaTA  for  CUBIC  SPLINE  method  ANO  GAUSS  QUADRATURE.  ARE  ALREADY  IN. 
rEADIS* 1U«)  NGaUS 
mGAUS1*NGaiiS+1 
m2*NGAUS/2 
nGAUS2sN2+  | 

PEAD (5>5GO)  (T(I),I*NGAUS2,NgAUS) 

PE AO (5, 560)  (M(l) ,1sNGAU32,MGAUS) 

DO  25  IQ»1 , Nil 
T(I«)s-TINgAUSI-IQ) 
as  miicdsmingausi-id) 

mRITLC6,56i)  (T(I),I=NGAUS2,NGAUS) 

WRITE (6, 562)  tw(l)  ,I=NGAUS2, NGAUS) 

5h0  pORMAl  (4E2().lO) 

5*1  fOHM A T (1X,aIII)sa*10(F10.3*1X)) 

5*2  pOHHA T(1X,amII)s**10(F10.S,1X)) 

L Rp  AD.  I N OAT  A*w*  ******  ■*****•***  ••A***'****************'*-**  A**  *************  ******** 

L MMaNMMHER  UF  CONTROL  POINTS  ON  THC  BLAOE  RADIAL  LOCATIONS  TO  BE  EVALUATED. 

C (hUSt  BE  AN  ODD  NUMBER) 

L M*Xlr«  MAXIMUM  NU.  OF  ITERATIONS  AT  WHICH  THE  ITERATION  IS  STOPPED. 

C XhmH'iB  RADIUS/T  IP  RAOIUS  RATIO. 

C X.1JB3PEED  COEFF ICIENT  (*VA/ (n*D)  ) . 

L 2a*N||MBER  UF  BLADES. 

C 3tGV*CAVITaTI0N  NUMBER  BASED  ON  THE  SHIP  SPEED  VA. 

c xx(It*normalizeo  radial  points  m bc  evaluated. 

L Bp TAC( l IsGEOMEtRIL  blade  SETTING  ANGLES  IN  DEGREES. 

t A(_FE  f I ) s ASSUMED  FFrFCTl  VE  FLOW  INCIDENCE  ANGLES. 

t Mp  Ib  no.  uf  terms  in  rnuRicR  sine  series  fur  gm.nf  to  he  greater  than  mf. 
c Np  ip  no.  uf  terms  in  fouricr  cosinc  stRiEis  for  in. 

c NpiLrN  is  nuibeR  nr  terms  (cven  no.)  for  filon  formaul*. 

C V*3HTP  IS  an  AdVANcF  speed  in  eeet/scc. 
c Ota  t3  propeller  diameter  in  inch. 

L InLOaO  FOR  NO  0*1*  FOR  VIVA(l),  JMl.D.Nt.V  FOR  USING  PREVIOUS  DATA. 

C IwRIt*  13  no.  OF  LAST  TTCWATIUNS  POM  which  THF  results  ARF  PRINTED. 

C Ep  Ip  INCREMENT  ratio  FUR  partial  OCR  T V A I t VES  FO  D ( CL ) /n ( ALFE ) . 


Preceding  page  blank 
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C GrGI  ( I * J»  K ) IS  NunnALIZFD  CIRCULATION  A3  FUNCTIONS  0|>  XX, SIGMA.  AND  ALFA 
t — — (t*  GAMNA/(VC*0>). 

u GrG»rGGI/(2.*PaI*R*VA) 

96*0(5,101)  fP 

96*0(5. ton)  MAxIT,HM,MF,NF,NFILON,IOLD,IwRITC 
96*0(5.1  G\>  OIa.VASHIP 
96*0(5,  IM)  XJJ,ZZ,3IGVA,XH,XXM 
mm1*hm+i 
mM2sMM-t2 
mmX=mm*2 

96*0(5, I U 1 > IXX(I),I3^,HM1) 

96*0(5, 10 1 ) (BtTAC(t).I=2,MMi) 
pE*0(5, 101 ) (AlFE(I),I*2,MM1) 

96*0(5,  1<M)  VS0LI(i),I=2,HHl) 
tFHULO.EO.O)  CO  10  722 
96*0(5. IU 1 j IVIVA(I),IB2,NH1) 

722  CONTINUE 
InO  FflRNA  T (81  to) 
lol  FORMAT (8F  ttf. 5) 

C Lrsr  THE  RLAO-IN  IUT*. 

wHI  TC  (<,,  1 80}  M*XIT,MH,MF,NF,IOLO 
wR!Tt(G,7l5)  01* 
wRITt(G,7l„J  V*SHIP 
wRITL(0,339)  NflLUN 
wRITLIb,  i 1<j)  ZZ 
wRI TC (b, I 85J  XH 
hRITC(0,U|)  XJJ 
wRITC(<*.  I 13)  StGVA 
9RITC(G,lt3>  (XX(IJ,TS2,MMI) 
wRITt(G.U/|J  »BFTAG(I),I*2,mnI) 
wRITL  ((>,  1 15)  (ACFE(I),Is2,MNi) 
wRI  TC  (•>»  1 77)  (S0LI(I)»Is2,HMi) 

TFllOLO.Ea.O)  RO  TO  724 
wRITC(0,723)  (VIVA(I),I=2,MM|) 

734  rONIINUE 

110  pORM*T(l3x, fNUNHCR  OF  BLA0ESs*,F3.0) 

111  pOHMAK 17X,*SPCFD  C0CFF.s*,F8.4) 

1 1 2 FORMAT (3X,*CAVITAT1DN  NO.  BASED  ON  VAs*,FH.5) 

1(3  FORMAT (5X,*XX (!)■*, 8(F10. 5, IX)) 

1 1 4 pOHMAT (2X,*BCTaG( 1)39,8(610.3, IX)) 

1 1 5 fORMAT(3X,*ALFc(I)»*,8(F10.5,1X)) 

177  fORM*T(3X,*SULI ( I )3*,0(F 10.5,1X11 
1b5  FORMAT ( 3X,*HUB/TIPs*,F6.3) 

1r8  FCIHMAT(tHt,l<IX,*MAXlT**,13,lX,»MMs*,I2,lX,»MF=*,J2,lX,*NF3*,I2,1X, 
X»IUL0**.I2) 

3^9  pORMAT  ( 10X,  »NF  IL(JNs*r  13) 

715  FORMAT (5X,*PR0PFLLER  0IAMET6R  IN  (NCH«*,F10.5) 

7 16  fOMMATISX^ADVANCE  SPEED  IN  FEET/SEC.**, FIO. 5) 

72S  FORMAT (3X ,*VIVa(I )s*. 8 (F 1 O.b.lX)) 


L This  13  A RUUTjNf  for  FINDING  CIRCULATION  FRUM  2-0  CASCADE  DAT*.#*** 
t.  IsIOsNUMBER  UF  OATA  points  TN  one  CL-SIGMA  CURVE (LESS  THAN  |«). 

C 3rGMTN  IS  MIN.  SIGMA  AT  WHICH  ALFAs6,  IN  THIS  CASE  IS  CALCULATED. 

L MANGtEsNUMUCK  0*  INCIDENCE  ANGLES  AT  ONE  BLADE  SECTION(LESS  THAN  S). 
t THjrSE  NUMBERS  MUST  BC  THE  SAME  FOR  ALL  XX (I). 
nO  120  IXXS2,MM1 

pE AD <5. 102)  XXX, MANPLt , IS IG, S I GM IH(  I XX  ) 
lfl2  fOHM A T(FIOa3,2ll0,Ft0.5) 

wRITCU,,!  t(,J  xxx, mangle,  1SIG 
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mPITC(G,665>  SlGMIN(TXX) 

1(6  fOHMAT(/,5x»*X«*,F5.3,5X,*NU.  OF  INCIDENCE  ANGLES  =*, 12,5X,*NO.  OF 
X DATA  POINT*  ON  ONE  CL  CURVE**, 12) 

6«S  fORMAT(5X,*MIM.  SIGMA  FOR  USE  OF  THE  MAX.  ALFAIlN  THIS  CASE  6 OEGR 
XpES)  **,F 7.4) 

L 3tGO(I,J.K)s  DATA  POINTS  FUR  SIGMA  WHERE  CL  DATA  ARC  READ  IN. 

C IN  ORDER  FROM  SMALL  SIGO  TO  LARGE  ONES, 

t ANGlfI)=INCIUENCE  ANGLES  USED  IN  2-0. 

L C(_  1 (7, J»K1=2D  LIFT  COEFFICIENTS  NORMAL  TO  VI. 

C Cf)l(t,J,K)=2l>  DRAG  COEFFISIENTS  PARALLEL  TO  VI. 
r*fl  120  IAMg*1,mANCLE 
pEAD(5»(0()  ANgKIANG) 

READ (5. tun  ISIGD(IXX»IANG,I), lat.ISIG) 
rEADIS#  tUl)  (CLI(IXX,IANG,I),Ia(,  ISIG) 

rEAD(5«  tUt)  (Lot (IXX.IANG. I), 1*1, ISIG) 

READ  (5,  lOt)  CGGGKIXXr  IANG.  D.Isl.ISIG) 

wRITC(6,11g)  ANGKIANG) 

wHITC(b,ll7)  (SlGDttXXrlANG.I)* 1=1, ISIG) 

WRITE (6,1  to)  (CLl(IXX,iANG,I),l3t,ISIG) 

WRITE (6,173l  (c01(IXX,IANG,I),lal,ISIGl 
WRITE (6, 8/19)  (GGGIdXX,  TANG, 11,1  = 1,  ISIG) 

8r0  FORMAT (IX, *GGGl (I,J,K)**,6F|0.5) 

1(7  FOKMAT(lX,*5IGD(I,J,X)a*,8F10.5) 

1(8  pORHATUX,*  CLl(I,J,K)a*,8F10.5) 

1(9  FORMAT (5X , *ANG la*, F5.2) 

tr3  fORMAUIX,*  CD((I,J,K)a*,8F10.5) 

120  CONTINUE 

L RpAO  IN  THE  DATA  FoR  CAVITY  THICKNESS. 
nO  360  IXXs2,HHl 
tFUXX.EU.2)  MRITE(6, 360) 
rEAD(5, tU2)  XXX,MANG1,ISIG1 
on  360  IANc.al,MANCl 
pE*D(5»  tUt ) ANG2(IANG) 

READ (5, tut)  (S:GC(IXX,IANG,I),Ia(,ISIGI) 
rE*D(5, tUt)(EE(IXX,IANG,I),I=l,13IGt) 
wRITE(6,36i)  ANG2UANG) 

WRITE (6, 362)  (SIGE(lXX,IANG,r),I=l,ISIGl) 

3*0  wRITE (6,363)  ( CF ( IXX , I ANG, I),lal,isiGt) 

361  FORMAT (5X, *ANG2**,F5.2) 

362  FORMAT (tX,*SIGc(I-J,K)a*,8Fl0.5) 

363  FORMAT (3X, *Ft ( I, J,K)a*,8Fl 0.5) 

360  FORMAT!//,  (X,*— — OATA  TOR  CAVITAION  THICKNESS——*) 

L CllRVp  FITTING  r)F  LAVITY  THICKNESS  BY  POLYNOMIALS. 

CALL  FI TC A V ( MM | , S IGF , EE , MANG I,ISIG1»ACC) 
oO  731  1=2, MM1 
nO  731  J=l,<i 

731  uRITE(6,73g)  ( ACC ( I , J,K ) , Ka( , 3) 

73O  fORMAT(1UX,*ACc(I»J,K)**,10F|0.5) 

C 2-D  rtA I A READ  IN  ARE  PASSED  ON  TQ  FUNCTION  CLCD ( I MM, SIGV , ANGLE, 1CLCD I . 


C CALCULATE  BClA«ATAN(VA/UMtUA*SMALLR)=ATAN(XJJ/(PAI*XX(I))) 
6 C»LC||LATC  ALF I»6t  T AG-ALFC-BETA . 

c calculate  ggam-ro-gfiag. 

oO  123  K 1 =2, MM ( 

C,GAM(Kl)«9g.-HcfAG(Kt ) 
pSl=XJJ/(PAl*XX(Kt)) 

RETaATAN(Psl) 

rETA(K1)*HeT*L0NVI 

*LF6(K 1 )»HE  TA6(K1J-MCTA(k1) 


c 
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ta5  »LFIC«n»fleTA6(Kl)-BeTA(Kl)-*LFC(Kl) 
wHITC (0*1 5bJ  (BETA(I)»I*C*HMl) 
i<;6  pt'WBAT(/»3X»'»«ETA(i)s*»«(F10.5»lX)) 
WBITCC0.157)  (ALF I ( I ) » 1*2*  MM t ) 
tsT  FORMAT t 3X,*ALFKl)a*,8(Fl0.5,lX)) 


C C*LCi|LATC  S1GV  FROM  AN  ASSUMPTION  FOR  THE  THE  FIRST  ITERATION.**************** 
C V | V A { i ) FON  ITeRA  NU.  greater  than  and  EQUAL  to  2 NAS  calculated  * 

t AFTER  EACH  ITERATION.  * 

t Vr V A f I ) FOR  ITeRA«i  IS  CALCULATED  BY  AN  ASSUMPTION  OF  VI*V*COS( ALF I) . * 

OO  J 26  ISs2,MM1  * 

XALF(IS)=CONV*ALFHIS) 

XPET(IS)sC0NV*BFTA(tS) 

VlVACIS)acoSlXALFtI3))/SIN(XBET(IS)) 

yVIal./VIVAdS) 

vvI2*vvi**2 

t2<>  SIGV(IS)»SIGVA*WI2 

wRITE(0,12b)  CSIGV(T)* I*2»  MM l ) 

128  FORMAT (3X,*SlGV(I)a*,H(FtO. 5, 1XJ) 


C FtND  GAMMAt  GGgUT)  DI STH IBUT ION .( NORMAL IZtD  TO  GGG/(0*VD).  * 

C FOR  THE  FIRST  ITERATION. 

HANsMANGLE 

manasman-i 

mANBsMAN-2 

no  130  NG=2,MM1  • 

AFSALFCINGj  * 

SC»SIGV(NG) 
wHITC(0,3«,5)  sc 

L — „-CA  VITA  r ION  NUMBER  CORRECTIONS . 

CALL  CAVlNn(NG,SC,AE»*CC.ANG2»ISIGl»SCN,EE,MANGl) 

wRITL(0.36«,)  SCN 

SR=SC/SC% 

SR=SQRT (SR) 

SC*SCN 

lCLCD*0 

CLS(NG)sCLcO(NG*SC,AC»ICLCD) 

iCLCDsl 

rDS(NG)sCLcUtNG*SC,ACf ICLCD) 

L C SUpRlIUTINE  CL CD  ALSO  INTERPOLATES  NORMALISED  GAMMA (CIRCULATION) . 

ICLCD*2 

130  GGG(NG)sCLC0lNG.SCrAC»ICLCD)*XX(NG)*VIVA(NG)*SR/Z2 
WRITL(0, 193>  (CLS( lit Is2»MMl ) 

WRITL(G,J9/)J  (CDS(I),I=2*MM1) 
hRITL(<*,131  ) (GGG(I),I=2.MMI) 

131  FORHAT(2X,*GAMMA(I)=*,10(F8.R,lXn 
228  CONTINUE 

RGG(D«0. 

nGG(MM2)an, 

t-  Find  COEFFICIENTS  OF  FOURIER  SINE  TRANSFORMATION  FOR  gamma  distribution. 

C PhKt)*IRANSF(irmEP  coordinatfs  for  xxci). 
pMI( i )=0. 
pHI(MM2)«PAl 
XALSJ.-XH 
*Hl=».5*(t.*XH) 
yHA*0.5*XA| 
nO  1 32  IPh32*mmJ 
*L*(XH1-AX(IPH))/XMA 
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r-  r-  r.  r r 


1*2  pHKXPH)»*tOSt*L) 

U$E  c^HIC  SPLlMf  METHOD  IN  COMBINATION  KITH  FILON  INTEGRATION  FORMULA**** 


FOR  integration.  * 

wFIHl«NFlLnN/2*l  * 
SPACL»PAI/nFILON  * 
nF11«NFILON*1  ' * 
OO  134  ICH* 1 , Nr  I i « 


134  9(ICH)*SPAcF*ClCH-l> 

Cl»*.E-7 

CALL  SPLINE^*!  ,GGG.0Y,S2,  S3.  S,SOUT,  331,332,  MM2,  NF  1 1 ,C1) 

ilO  145  ISFK*l»MM 

aRG>ISCK*SPACF 

aMH*FALFIAR6) 

bMH*FBL T l ARG ) 

RHH*FGAM(ARGJ 

AFFsANh«(SoUT(l)-Sn(IT(NFIl)*C03(ISLK*PAI) ) 

92K=0. 

j2KA*0. 

r»n  510  K2»  ( , NF  IH1 

t2A»K2*2 

t2=K2*2-1 

S2KsS2K*Sni|T(l2)*3IN(ISEK*3(l2)) 

IFIK2.UU. NFIH1)  GO  Til  510 
s2KAsS2KA*SOUT(I2A)*SIN(ISEK*S(I2A)) 

5(0  CONTINUE 

p32K38HH*S2K 

f;S2KA=GMM*s2KA 

RMllSEK)=SpACr*(AFF+B32K+G32KA) 

l«6  CONTINUE 

rMII3EK)=  gMCISFK)*?./PAI 
WRITE (6,  I/I7)  ISFK,GM(ISEK) 
ln7  F0RMAT(5X,3HGM(,I2,2H)3,E14.7) 

225  CONTINUE 
1«5  CONTINUE 


Calculate  induction  f actons. ********************************************* 
calculate  small  ia  and  in.  * 

THEN  CALCULATE  INA  ANO  JNT.  * 
no  150  id*2, mhj  « 
XBE=(0LTAG(ID)-ALFE(IU))*CONV  * 
TBIs I ./T  An ( XBF)  * 
TBI2»rBI**2 


E PhI(|)— *X(I)  HAS  BCEEN  DONE  BEFORC. 
wFl>NF*t 
r»n  151  nx*  ( , nf  1 

NX A*NX“l 

tFINX.GE.2)  GO  to  3P9 
nO  326  NI*(,NFll 

t Change  of  vahi abuts  from  sm  in  xxpii). 

XXP»XHl-XHA*COStSfNT)) 

xiNo*xxp/xx(in) 

XlI*l./XlNf) 

»B1«ZI*(XinD-1.) 

aB2*AB1*TBI*XIi 

e Calculate  ai  ano  as,  then  hi  and  02. 
vpp*xii*thi 
yPP2*yPP**2 
qC*l.+TPP2 

pO»l.trHl2 
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nC3*90HT(Rc) 

hDS«SQMT(Bd) 

nKl*bC3-U0s 

oK2*lB03-t.)*(BCStl.) 

,>K3»(B03tt.)*(BC3-l.) 

*DK*a.S»ALnGl0K2/0K3) 

*A1«OKI-A0k 

4A2a-AA! 

oFl*8D/BC 

nF10«0Fl**u.25 

oF2sTBX2/Ro**1.5 

pGAl*EXP(ZZ*AAl) 

FGA2»eXP(Zz*AA2J 

EH»l./lEr.Al-l.) 

E2I*l«/lfpA2“j.) 

TFIXXNO.lt. 1.)  GO  TO  327 
AEG2*AL0G( i .♦  C2I) 

«B2»0F10*(E2t»0.!I*0F2*AEG2/ZZ) 

XlACiDfNfJs  AB2*(1.+BS2) 
xTTCXIJfNT)  - AS  1*682 
rO  TO  320 

3g7  *EG1  *ALl)(»(  j ,♦  EU) 

rB1«0F10*(eH-0.5A0F2*AEGI/ZZ) 

*:a(xo,nt)3-ab2abbi 

xir(XO,Nn»-ABlA(l.+BBl) 

326  c'1NTXNUE 
32‘)  rONTINUE 

C2KA=0. 

c2XTao. 

CZKAAsu. 

C2KAT»0. 

f)0  350  TCHal  r NrIH  1 

t2»ICB*2-1 

t2A*XCB»*2 

cMT|aC0SlNXAAS(I2)) 

C.’KAsC2KA*xIA(iO,X2)*CNTl 
c2KT3C2KT  + xnU0,X2)*CNTl 
tFUCB.EU.nFIHD  CO  TO  350 
rNT2«C0SlNxA*3(I2A)) 
c2KAAaC2XAA*XXA(I0,I2A)*CNT2 
r2KAT*C2KATFXIT(I0»I2A)*CNT2 
3s 0 C(INTXNUE 

C2KAsC2KA-.5*(XlAf Xn,NFIJ)*CnS(NXA*PAn+XIA(IO, 1)) 
c2KTsC2K!-.5*(XlT(Xn#NFin*C0S(NXA*PAl)*XITU0,n) 

AR(iaNXA*SPACE 
BNHrFBtTURGJ 
(•NHaFGAHl  ARC) 

SlNA(I0*NX) aSPACC* (BNH*C2K AtGNH*C2K A A ) /PAX 
SXNT UU# NX) aSP ACE* (BMH*C2KftGNH*C2KAf)/PAI 
rFlNX.GE.P)  3lNA(J0,NX)a  2.*SINA t 10, NX 1 
lSi  tF(NX.GE.2)  3InT(X0»NX)s  2.*SXN T ( ID, NX 1 
ISO  fONTXNUE 

t FtNO  HHA(I)  and  HHT(I),*»**»*************«******»*»**«**««»»«*»««*»*****' 
oil  169  MAS2»NMJ 
jPaSIN(PHlcMA)) 
oO  160  MBs j, MM 
pHaMB*PHl ( M* ) 

SPM3SIN(PM) 

fPMSCIlStPM) 
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HHACHA»HB)aO. 

hHT(HA.HB)«0. 

no  161  HCa l » MBJ 

mC  A*HC“ I 

pH  I *MC  A*I*H  J ( HA  ) 

CPHl»C03iPHlJ 

HHA(HA.MB)sMHA(HA,HB)*SINA<HA.MC)*CPHl 
1*1  nHT(HA.Hb)sMHT(HA,MB)*SINT(MA,HC)*CPMI 
HHA  t HA,  HB ) aHHA (HA  * HH ) *SPN 
HHT (HA  # HB) aHMT ( MA  » MB ) *SPH 
mB2*MBa2 

nt)  162  HPa*B2,NFl 
"PAaMP-t 

<|NP*SIN(MPa«PHI(HA)) 

CPS*CPH*3Hp 

HHA(HA.HB5«HMA(HA,HB)+CPS*3INA(HA,HP) 
1*2  MHT(MA.HB)aHHT(MA,MB)tCPS*3rNT(MA,MP) 
pSPaPAI/Sp 

hMACMA,MB)*PSP*HMA(MA,MR) 

HH  r ( HA»  HB ) sPSPaHMT ( M A * MB ) 
lfcO  CONTINUE 
1&9  cONJINUE 


t.  NpHTnN  ITEHAlIoN  LOOPa******************************************************** 
L RfADy  TO  EXENClSE  NEWTON  METHOD. 

t RfNHjTE  ALFEtl)  Ano  GH(t)  BY  A N£N  SET  OF  VARIABLES,  * 

C XN(I)»ALFE(D,  i a 1 , HH  AND  XN  ( I ) aGH  ( I J , IaHHl  , HHX  . • 

fTERAai 

TTERAHaHAXiT-lwRITC 
RB7  CONTINUE 

n,'ITCCG,66l)  ITFPA 

6fcl  fORNAT(///,  10X,*—— ITERATION  NUHBER— — A,  13) 

HO  525  INEal , Hh 
fNElalNE+I 

aLFEO( INEl )*ALFE( INEt ) 

RMU(INE)»GH(tNf) 

SlGVO(INEl)aSlGV(INEI  ) 

Vi V AO( INEl )SVIVA( INEl) 

535  *N( ine ) »alff IINFI ) 
r>n  526  INEsHHl,HHX 
TNEPalNE-HH 

5j6  yN ( I HE ) aGH ( INfP) 

wRITt(6,i»3/|)  (XN(I),t  = l,HHX) 

S3O  fORHAT (tX,#XN(I)3«,12(F<>.5,lX)) ’ 
r»0  170  LUaj.HNi 
HAVA (LU)»0# 
hTVA(LQ)»0. 

00  171  LPa 1 ,HM 
XPAG*LP*SH(LH) 

HAVA(LU)3HAVA(L(l)tXPAG»HHA(La,LP) 
t IX  wTVA(LU)3NTVA(L0) ►XPAG»MMTtLU*LP) 

HAVA(LU)»HaVA(L«)/XAL 
I7O  HTVA(LU)awTVA(L'J)/XAL 

TFHTERA.LT. ITfRAN)  GU  TO  22 h 
HRITE(6,25/|)  (MAVA(I),Is2,HH() 
mR  I TC  <6*25s>  (WfVA(I)/I *2  * HH 1 ) 

2j6  CONTINUE 

2S0  FORHAT  UliX,*HAvA(I)a*,B(f  10.3,  IX)) 
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2*5  FORMAT (tux, **TV*(I)**#B (FI  0.5# ix)) 


r 


C F ( I ) 13  AN  array  or  FUNCTION  AND  P(1,J)  IS  PARTIAL  DERIVATIVES. ••••*«**•• 
00  530  IFNalrMN  * 

tFNIaIFNft  • 

RCA(lFN)s(Br TA6( IFNI )-ALFE( IFNI)) ACUNV  * 

?8G ( IFN) *t  an i bra ( IFN) ) • 

pXJ(IFN)»PaI*XX(IFN1)/XJ^ 

530  f(IFN)3TBG(IFN)*(PXJ(IFN  )-WTVA(IFNt))-(l.tMAVA(IFND) 
on  531  lFNaMMi,MMX 

r.MAso. 

TFN]sIFN-MMfl 

oU  532  IFMai,MM 
lG*IFM*PMI(IFNl) 

532  rMAsGMA^GM(IFM)  *SIN(AG) 

531  f(IFN)»GMa-GGG(IFN1) 

nRITE(U,533)  (F(I),I31,MMX) 

533  FORMAT (tX,*F(I)3*,12(F9.5#lX)) 

l computation  of  partial  scrivatives. ***********  ,***»*************«******** 


r»0  535  IPAsl#MN 

* 

oO  535  IPBal#MM 

* 

rFdPA.EU.lPB)  GO  TO  536 

* 

p(IPA#IPB)sO. 

* 

R0  TO  535 

536  TPA)sIPA*t 
CBGA»COS(BGAtIPA)) 

CRGA2aCBGA**2 

aXY  =-(PXJdPA)-*TVAdPAl))/CflGA2 

p(IPA.IP9)aAXY 
535  CONTINUE 

r»0  537  IPAal#MM 
TPAJrlPAtt 

r>0  537  ipbsmmi,mmx 

TPBSaIPR-MM 

SSL  a" f HG ( IP A ) *MMT(IPAt#IPBS)*IPSS/XAL 

537  p ( IPA » IPB) sSSL“HmA ( IP A t , (PBS) *IPBS/X AL 

t rpG(r)  and  pxj(i)  are  alrlaoy  calculated. 

oO  538  IPAaMMl , MMX 
rPASalPA-MM 
no  538  IPBalfMM 
tFIIPAS.EO.IPH)  GO  TO  538 
P(IPA,IPB)bO. 

538  CONTINUE 

t 0(CL)/D( ALFE)  NEEDS  A FINITE  DEFERENCE  METHOD. 

OO  539  ICLa2»MMl 
0lFaCP*ABS(ALFL(ICL)) 

*FXaALFE(ICL)»OlF 

AFYaALFE(IcU-0lF 

SC*SIGV(ICL) 

CALL  CAV1N0(ILL»SC#ACX, ACC, ANG2, ISIG1 ,SCNX,EC.MANG1 ) 

CALL  CAVlNOdCLrSC, ACY , ACC , ANG2, ISIG1 . SCNT ,EE #MANG1 ) 

jCXaSCNX 

jCYaSCNY 

SRXrSC/SCX 

jRYaSC/SCY 

gRXaSUHT (SRX) 

BRVaSOHT(SRY) 

TCLCD32 

fiCGXeCLCD(iCL,SCX.AFX#  ICLCD)-*SHX 
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G 


c 


* 


ft 


G 


rGGY*CLCU(ICL,SCY. AFY,ICLCD)*SRY 
TCLA»ICL-l 
TCLAM*iCL4*MM 
XVISXX(JCL)*VIV*(ICL)/ZZ 
5y>  p(ICLAH,ICL*)a-CGGGX-GGGY)«XVl/(2.*0IF) 
n0  540  IGsi.MM 
tgi»igm 
pHNaPHi IlGl> 
nU  540  IK* 1 r MM 
yKP*IK*PHN 
SXK*3IM(XKP) 

TGM*IGAMM 

TKH*XX-fMt1 

5fl0  p ( IGH  t IKH)  aSXK 
no  541  IHTalfMMX 
TFUTEWA.LT.lTfRAM)  GU  TO  582 
5/|l  HRlTC(0,5n2)  (P(IPT,t),I*l,MHX) 

5a2  pORMATUX,*PlIf  J)s*.  t2(F9.5»lXJ) 

5fl2  c«NT^NUE 

CALL  DCTERM(P,NMX,OETBO) 
nO  543  lPF*i/MMX 
nil  544  ipg*1»mmx 
n(lPG*IPF)*P(IPG,IPF) 

5afl  pIIPGf IPF)3-F(JPG) 

CALL  DETEHH(P*mMX»OETE(IPF)) 
no  545  IPG* I. MMX 
5a5  p(IPG,IPF)*ttClPG,IPF) 

5a3  cnNT!NUE 

nO  546  LS*j,mmx 

5a6  *NlLS)»XN(L3)+orTE(L3)/OETBO 
tFUTEHA.LT. UE«AH)  GO  TO  580 
wPITC.16,548)  CXN(I),  1=1, MMX) 

5a8  F0HMATCtX,*XN(i)**,l2(F4.5f IX)) 

5«0  CONTINUE 

00  547  L A*2»  MM I 
LAA*LA-t 

ALFECLA)=Xn(LAA) 

1 A AM»rL A A + MM 

547  BM»LAA)*XN(LAAM) 

C C»LC||LATC  UCTAI#ALFI,ALFC»VIVA  ANO  SIGV.**A 
C NfN  pAVA  ANO  HTVA— — - — 
r»0  550  MAT*2»MMl 
W A V A ( MAT ) *0 . 

MTVA(MAT)3>J. 
on  551  MAS* 1 * Mm 

xpag*mas*gm ( ma$ ) 

pAVA( MAT) *KA V A (MAT ) ♦XPAG*HNA(MAT»MA3) 
551  mTVA(MATJ*NTYA (Ha  T) tXPAGAMMT  JMAt»MAS) 


WAVA(MATJ*MAVA(MAT)/XAL 
5*0  MTVAtMAT)swTVA(MAT)/XAL 

nn  180  KF*2» MM1  4 

XK1SPAI*XX(KF)/XJJ  * 

PETB*C 1.+M*VA(kF))/CXK1-MTVA(KF ))  A 

qf TAI{KFJ*C0NVlAATAN(»ETH)  « 

4LFICKF)*8eTAI(KF)-8ETAfKF)  * 

*B«8CTA(KF)*CONV  « 

»BI*8ETA1(kF  }*C'JNV 
p8*CU3(XB)/SIN(XB) 

CHIsl./CU3<X0I) 


I 


t 
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vIV*CKF)»(cH-wtV*(kF))*C8I 
1«0  sIGVCKF)*S|0V*/VIVA(KF)**2 

TFIITERA.lt. ITl-RAW)  GO  TO  227 

mRITEIG,I8|)  3 

wRITCCG,196)  ( ALFGC l ),I*2,MMl) 

wRITCCG, 19j)  (BFTAJ a)#I*2*MMl) 

gS I TL (0» 1 89)  ( ALF I ( I ) , I 3£ , MM ( ) 

nRITCCC.183)  (VlVA(I),I*2,MKl) 

wRITCCG, Ifl/i)  (SIGV(I).I32,HH1) 

2g7  CONTINUE 

t«l  pORHATtlOX,*— — NEW  VALUES  OF  ALFE,  V I VA.ETC  .— — *) 

183  fORmaT(3X,*VIVA<I)s*,8(F10.5,1X)> 

108  FnHHAT(3X,*SlGV(I)3*,S(F10.5,lX}) 

1 «9  FORMAT (3X,aALFI(I)3*,8(F1 0.5, IX)) 
tq6  fORMAT13X,*ALFG<I)»*,8<F10.5,1X)) 

197  FORMAT (2X,*BCT  Al (I)s*,8(FI0.5,1X)) 

c calculation  of  propeller  fupces  and  efficiency. a************************* 


l Calculate  unif/va,  alfinf  ano  ggg.  * 

no  JRO  x r=2,MM j • 

AEaALFC(KT)  * 

SC»SIGV(KT)  • 


CALL  CAVlNo(KT,SC,AEr ACC,ANG2,ISIG1,SCN,EE,MANG1) 

jjRI (KT )sSC/SCN 

SRI(KT)3SQrTCSRI(KT)) 

tFUTERA.LT. ITERAR)  go  TO  775 

wRITCCG, 365»  SC 

wRITEtG,36b)  SCN 

3fc5  fORMAT (5X,»C AVI  CATION  NO.  BEFORE  CORRECTIONS*, FlO. 5) 

3*6  fORMAT(5X,*  CAVITATION  NO.  AFTER  CURHECT IONs*,F 1 0 .5) 

775  CONTINUE 
SC=SCN 

TCLCDsu  « 

CLS(KT)sCLcOlXT,SC,AC*ICLCD) 

tCLCDal 

C0SCKT)sCLc0tKT,3C.AC*ICLCD) 

tCLC0*2 

rGG ( K T ) sXX (K  T )*VI V A(KT) *CLC0(KT , 3C,AE,ICLCD)*3RI(KT)/ZZ 

I9O  CONTINUE 

tFUTERA.LT. ITERAR)  GO  TO  228 
mRITEIG, 1 90)  ( GGG ( I ) , 1*2, MM t ) 
wRITCCG, I93)  (CLS(I),I=2,MM1) 
wRITECG,190)  (ci)SU),I=2,MMJ) 

193  fORMAT(5X,*CL(I)s*,8(F10.5,1X)) 
l9'l  FORMAT (5X,*LD(I)s*,8(F10.5,1x)) 

1 q8  FORMAT (2X,*GGG( I )3*,8(F 10.5,  IX)) 

2?«  CONTINUE 

C OaTA  FOR  CUBIC  SPLINE  MFTMUO  ANO  GAUSS  QUADRATURE. 

XCT(1)«-1. 

XCT{MM2)s». 

fcttd»o. 

FCP(1)»0. 

FCT(MM2)*n. 

FCPTMM2)s,t# 

fTD(I)>0. 

FPD(1)sO. 

fTDTMM2)30. 

FPD(MM2)*0. 

bESsOIA*.5/12.*VA3IIIP»I.E5 

rE1»I.L6 
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nE2»l.E8 
on  310  ICT»2,MMl 
XCTaCT)*(xX(ICO-XH1  J/XHA 
xcoN«XHitxH**xc r ( ICT) 
nlEla  BEfAiC ItTloCUNV 

rBIE»CUSlBlf l) 

«BlE»SIN(H|rt) 

FCaVlVA(ICT>**2*SOLT(ICT)*XCON*XMA  *SRI(ICT)**2 

FCT(XCT)3FC*UBIE*CLS(ICT>-SBIE*CDS(ICT)) 

FCP(lCn»Fc*XtnN*(3BtE*CL3tTCT)tCBIEAC0SUCT)) 

C RpTNoCOS  NUMBER  BASED  ON  PROPELLER(RCS)  ****FRXCTION  DRAG** 

C RfYNdCDS  NUMBER  BASED  ON  PROPELLER  RADIUS  AND  ADVANCE  SPEED— 
t.  — RES«R»VASHIP/NiJ. 

C RpX  tS  A LOCAL  REYNOLDS  NUMBER. 

VlVA2*VXVA(XCT)**.? 

pSD«VIVA2*sOH(Icr)*XX(ICT)*.5/XZ 

9EXaVIVAtlcT)*SOLI(XCT)*2.*PAI/ZZ 

BEXsHEX*XX(ICT)*RES*SRI(lCT) 

TflHEX.LE.nri)  CFa0.664/SQRT(REX) 

TFlREX.GE.nr2)  Cf30.0i0/REX*«0. 142857 

CFa0.044/R£X**(j.  166666667 

pSDF3PSD*CF*XH» 

nETAGR»BETAGllcT)ACONV 

SBG»SIN(BETAGR) 

rBGsCOBtUFTAGR) 

FTD(lCn*Fc*CF*SBC 
3(0  fPOCICT)3Fc*XLON*CF«CBG 

CALL  SPLlNE(XLT»FCT,DY,32»S3,r,OFCT,SSl.SS2,MM2,Nr,AUS,Cl) 

CALL  SPLINE tXLT»FCP,DY, 32# S3, T,OFCPfSSl.S32. MM2, NGAUa.Cl) 

CALL  SPLINE t XLT • FTQ,DV  ,32,33,  T,UFTD,  SSI  ,SS2,  Mi2,  NGAUS-, Cl) 

CALL  SPLINf tXLT ,FPO, OY , 32, S3, T, OF PD, SSI, SS2, MM2, NGAU3#C 1 ) 

CT=0. 

rP«o. 

CTD«0. 

cPOaO. 

nO  311  ICP«1,NRAUS 
cT*CT+OFCT(ICP)*n(ICP) 

CPaCP+OFCP(XCP)AN(ICP) 

cTO*ClDtOFTOlICP)*N(TCP) 

3| 1 rPUsCPUtUFpOlIcP)*W(ICP) 

XLAMaX JJ/Pt I 
f T*2.*CT 

rP»2. *cp/xl*m 

FFF*CT/CP 

CTD*CT-2.*cTD 

CPDaCP*2.*cPD/XLAM 

pFFD»C  fD/CPO 

XJJ2*XJJ*«2 

xJJ3*XJJ2»xJJ 

XKT»CT«XJJ2*PAI/B. 

xPT«CP*XJJ3/16. 

XKT0*CrD*X.;J2»PAI/n. 

XPTD«CP0*X./J3/1G. 

Tf (XTERA.LT.tTERAN)  go  TO  5«5 
wRITE(0,3t2>  LT  , XKT,C TO, XK TO 
wRITLtO,3t3)  LP  ,XPT,CPD, XPTD 
»RITE(fa,Sl4)  EFf,EFFD 

3(2  FORMAT ( 1 $X ,* THRUST  TOtFF IC IEM Ta* , F 1 0.5, 5X , *K Ta* ,|-  J 0 .5, 3X. *C f WITH 
XfHICTION  DrAC**.F  t0.5»2*,*" TO»*,F 10.5) 

3|3  FORMAT  UBX,*P<JWFR  COCfFlCXENTa»,F  t0.5,5X,*"Pa*,F10.S,3X,*CP  NIIH  F 


1 


XqICTiON  DRaG3*,F10.5*2X,»KPOs*,F10.5) 

3(0  pORMA f C22X, ^EFFICIENCY  CT/CP=*,FI0.S,5X,*EFFICICNCY  WITH  FRICTION 
X0RAG**»FIn.5) 

5*5  CONTINUE 

TFUTEH*.Ea.M*XlT)  GO  TO  999 
tTERASITEPa+1 
nO  755  IOk-2,mm1 
tuka«iuk-i 

GMUUKA)3XxH*GMO{IOKA)t(l.-XXM)*GM(IOKA) 
»LFE(I«JK)sxXM*ALFF.O(IOK)tU.-XXM)«ALFECIOX) 
aLFI(IUK}sbFTAG(IOK  J-BETA (IOK)-ALFK(IUK) 

7s5  BETA  I C I OK ) -Ul'T  AtlOK)-ALFI(IUK) 

TFIITEHA.Le.ITLHAW)  GU  TO  230 
mRITEIG,1B2)  ( ALFE( l )»Is2»MHi) 

1r2  FORMAT (3X #*ALFf II)=*»H(F10,5.1X)) 

2j0  CONTINUE 
R0  TU  987 
999  a TOP 
pNO 

*» 


r 


FUNCTION  CLCUCI.S.b.ILD) 

C THIS  SUBROUTINE  INTFRPOLAIES  CL. CO  AND  CIRCULATION  FROM  INPUT  DATA  AT  X.  SIGMA 
C ANO  ALFA.  OCPCNOING  ON  ICLCt)*0,l,2. 

CONNON  PAI,  CONV,  CONVI.SIUMINCIO) 

COMMON  SIGOUO.S.IS).  CL1  (10.5.15)  . ANG1  (5)  »C01  (10.5.  IS) 

COMMON  MANGLE.  ISI6.GGGIU0.5.1S) 

DIMENSION  PC  3,10), S(  5.10). STUN CIO, 231. ALE2 (I0),AC10) 

C I SPECIFIES  XXCI). 

i S II  CAVITATION  NQ.  FOR  ALG2  TO  BE  EVALUATED. 

i A is  EFFECTIVE  flow  incidence  angle  FOR  ALG2  to  he  evaluated. 

CHISOSU.UI 
ISIGA*ISIG-1 
m l an ANGLE- t 
mansmanglf 

IFcS.LL.SIGmIN(I))  mansmi 

DO  | K31.MAN 

ITa3 

IFcS.LL.SlGDlI.A.m  GO  TO  10 
TFIS.GE.SIGOCI.K.ISIG))  GO  TO  11 
IFIS.LL.SIGOU,*,?})  GO  TO  10 

ifis.gl.sigou.a.isiga))  go  to  11 

21  033IGD(l,f«,in-S 
IFlO.GL.Ii.)  GO  Tu  20 
IT*If>l 

IFllT.GE.I3IG)  STOP 
GO  TO  21 
20  13c*lT-2 
isb»it-i 

ITTalT 

ISlalT+1 

22  PCI.) )«3IGDCI,K. ISC) 

Ptl,2)*SlGD(t,K»ISB) 

Pf 1.3)sSIG0Cl,K,ITT) 

P(t,«)»SlG0Cl,K.I3l) 

IF(ILD.EG.i)  GO  TU  30 
IF  ( ILD.EU.2)  GO  T'J  32 
P(2,i)=CL1 CI.R.ISC) 

P(2»2)*CL1 (I,k , IS8) 

P(2,5)*CL1CI,k,ITT) 

PC2. A )®CL1 (1.K.1S1) 

GO  TO  31 

30  PC2.1 )3CUt (l.K.ISC) 

P(2,2)*CUt CI#R,IS8) 

P(2,3)=CDtCI,«,lTT) 

P(2, A)*CUt (1.K.1S1) 

GO  TO  51 

32  PC2.1 )=GGGIC l.K.ISC) 

P(2»2)*GGGKI.k,IS0) 

P(2»3)*GG6I (I.K.IT1) 

P(2,«)»GG6Ilt,K»ISl) 

31  CONTINUE 

DO  55  KS*l,4 
55  P(3,KS)«V.902 
Naa 

M«3 

CALL  LSGUARCP.N.M.A.CHISD.STnR) 

ALE2CM )■*( 1) 

MM*aM-l 

->0  5<*  KMSt.MMA 
KMlaMM+t 


1 

m 


H 


A 
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I 


56  *LE2<kj,a'-£21K)+6(Kki1J*5**KH 
SO  TO  1 

10  ISC* l 
130*2 
ITT*3 
isi*a 
CP  TO  22 
U ISC*XSiS"3 
ISB*XSIC-2 
ITT*lS16-t 
IS1*ISIK 
GO  TO  22 
1 CONTINUE 

C NO*  FOUND  ALE21NANGLE)  CuRRESPONOING  TO  ANGi(*ANGL CK 
u THUS  OETERHINP  LLCD  AND  CIRCULATION  CORRESPONDING  TO  A. 
DO  57  KL*  t >Z 
DP  57  KN*1 > NAN 

IFCKL.tQ-l)  aiKL,KN)3ANGllKN) 

IF(KL.C0.23  aiKL/KN)3ALE2tKN) 

57  IF(KL.EB.3)  utKL,KN)30.002 

*ANA*HAN-l 

nanB*han-2 

CALL  LSQUAR(Cl,HAN.HANA,A.CHI3Q»STOiO 
CLCD*A(1J 
DO  58  KT*1,HAN6 
K T 1 **  T+ 1 ' 

58  CLCD*CLCU+ACMlJ*a**FT 
RETURN 

END 

«« 
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SUBROUTINE  FXTC*V(MM1 fSE,EE,Ml,n#ACL) 

nIHENSlON  sfU0#5,5),EE(l  0,5,5),  ACC ( 1 0 . 5,5 ) ,P (50 , S) , 9 (50 , 5 ) 
r>0  1 1=2, H*1 
r>0  I 

nO  2 K»t#U 
nO  2 L=I»Il 

tFIL.EQ.1)  P(K,L)=EE(I,J,K) 

TFIL.EU.2)  P(K,L)=1. 

2 jFIL.EO.3)  P(k,L)*-SC(I»J»K) 

<-ALL  DET£»m(P,I1,DB) 

nn  3 IU3J,I1 
no  « LP*1 » 1 1 
0(LP. IU)=P(LP,I0) 

1 p(LP,ID)=Se(I,J,LP)*EE(I,J,LP> 

CALL  DETERm(P,II,DC) 

*CCCI.J»ID)=UL/0B 
nil  5 LPs1,TI 

3 p(LP,IO)3Q(LP,ID) 

3 CONTINUE 

I CONTINUE 
RETURN 


* 


SUBROUTINE  C A V I NO  ( NG,  3C  , AE  , ACC  , A2,  1 1 , SCN  , EE » H 1 ) 

DIMENSION  *CB(5) 

DIMENSION  *CC(  lG,S,5)rA2 (10), EE (10,5,51, XX (5) , p (50 ,5) , Q (50 ,5) 
pRl=l.E-3 

TFUsO 

mIUHI-1 

2 tFUE.LE.A2(U)  GO  TO  1 
TFU.GE.M1)  GO  TO  1 

T*I  + 1 

rO  TO  2 

1 continue 

K=t-2 

TFU.E0.J.n*».I.EQ.2)  K = 1 
TFU.EO.MlA.on.I.EO.Ml)  K=Ml-2 
k3sk>2 
t FjNO  EEX. 

„0  3 KI»k,k3 

L * MUgT  BE  LESS  THAN  ONE. 

a*;  5 

jO  qsA-O.OOl 

TFUFL.EU.I)  NRITE(6,O0) 

00  pL’RM  A T (IX, — C BECAME  NEGATIVE-—-*) 

TFUFL.EU.j)  STOP 

CALL  FANC(pF,B,SC,ACC,NG,KI, Jt) 

IFIPF.GE.O.)  GO  TO  11 
4=8 

rO  TO  10 

[1  pfWlNUE 

aAR=A 

A=« 

rsAAR 

CALL  MOSEC(A,U,FRl,ER2,X,J,NG,KI,SC.ACC,Ii) 

t Convert  c(OR  x)  into  sigma. 

pSl=2.-X 

pS2=2.-2.*X 

pS3=PSl/PS2 

pS«=PS3**2 

gIGMAsSC*P34 

3 *XlKI)=SlGMA 

C Wp  H*VE  CAVITATION  NO.  AS  A FUNC.  OF  ANG2. 

C PyT  TN  A POLVNoMlAL  FORM. 
nO  20  13=1,3 
no  20  IT*1,3 

20  p(IS,IT)=A2(K*lS-l)**(IT-l) 

CALL  0ETERm(P*3#DH) 

nO  21  10*1,3 
nO  22  IE*t ,3 
r(IE»ID)=P(IE,IO) 

2 2 p(IE»IO)=XX(K  ♦IC-l) 

CALL  0CTERm{P,3,0C) 
aCH( IO)*OC/OB 

nO  23  I E=l,3 

23  p(IE.IU)=G(ll»ID) 

21  CONTINUE 
gCN=ACB(l) 
nO  20  11=1,2 
TLl=IL+t 

20  gCH*SCN»ACR( IL I )*AE»*IL 
(?F  turn 
fnd 
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SUBROUTINE  FANC(F.CE'SC»ACC'NG,KI,I1} 

rtlNENSION  aCCU0,5,5) 

QH31aACCtNG.KI,n*CE»ACC(NG,Kl»2) 

QMS2aEL+ACc(NGfKI,3) 

pHSaHHSt/HH92 

»A1«2.-EE 

»A2a2.-2.*er 

A A3>A A 1 / A A2 

4 AR*AA3**2 

XLHS*SC*AA/| 

paXLHS-RMS 

ps-F 

RETURN 

gNU 


SUBROUTINE  MOStCCA,B,ERl,ER2,X,J,NG.KI,SC,ACC#Il) 
DIMENSION  *CC(10,5,5) 

J*0 

Xt«A 

X2*» 

A J»J>1 

rFU.GC.8ii,,)  GO  TU  8 

CAUL  FANC(pFXl,Xl,SC,ACC,NG,Kl#Il) 

CALL  FANC(pFX2,X2,SC#ACC,NG,Kl,Il) 
x3*Xl+(X2-xl)*PFXl/(PFXl-PFX2) 

CALL  FANC(pFX3,X3»SC,ACC,NG.Kl,Il) 
jF IPFX3) 1 ,2*3 

1 X2MX3 

X I = « 1 

tFlA-H)lU,lO»U 
,0  y = X j-ER I 

TFIY.LE.V.)  ra«. 
rO  TU  12 
ll  v*X3*ERt 

t2  c*LL  FANC(pFY,Y*SC,AC:*NG.KI,Il) 
tFlPFY)  5,2,2 
3 xl=X3 
X2»*2 

rFlA-U)  20,20,21 
2<>  78X3+EN1 
rO  TU  22 

21  2*X 3-EHl 

22  CALL  FANC(PFI,Z»SC,ACC,NG,KI,I1) 
rF  t PFZ )2,2, 5 

5 rO  TO  'I 

2 pPa  ABS(PFx3) 
tFIPP-CR2)  6, 6, A 

6 x = X3 

rO  TO  7 

8 wRITL(G,9)  J 

9 fORMAT(1X,2HJ*,I3) 

STUP 

7 RETURN 
pNO 


E 


r 


L- 


C> 


ft 


ft 


o 


I 


subroutine  detfrm  ca.n.dj 
L DETERM  REVISED  02-28-73 
REAL  N 

OtMCNSIUN  A(.50.50),5AyCA<50,50) 

IF  CN  .EO.  I ) GU  TO  46 
Cal. 

NM  * N 

1)0  ’ j * l»NN 
00  *>  1 * l»NN 
9 S4VCACI. J)  * A(I#J) 

K a 1 
GO  TO  13 

12  « a K ♦ t 

13  I 3 4 * t 

L * K 

GO  TO  17 

16  I 3 I ♦ 1 

17  IF  CABSCSAVEAII.K))  .ST.  A8S ( SAVEA (L.K ) ) ) L 3 I 
IF  (I  .N£.  "N)GO  TO 

IF  CL  .EG.  K)GU  TO  28 

J a K 

C R0«  INTERCHANGE 

GO  TO  23 

22  J a J ♦'! 

23  SAVEKJ  = SAVEa(k.J) 

SaVEA(K.J)  3 3*VEACL»J) 

Sa'/EACL.J)  » SAVEKJ 

IF  (J  .N£.  nN)GO  TO  22 

C 3 -C 

28  I 3 K 1 
GO  TO  31 
30  I 3 I ♦ t 
3l  CONTINUE 

IF  ( SA VEA (K « K ) .£0.  0.)  GO  TO  48 
M 3 SAVEA(I.K)  / SAVEA CK.K) 

SAVEACI#K)  3 0. 

J 3 K ♦ t 

GO  TO  3o 

35  J s J + t 

36  SAVEACIfJ)  3 SAVEA(IrJ)  - M * SAVEA (K.  J) 

IF  (J  .NC.  NNJCO  TO  35 

IF  Cl  .NE.  nn ) GO  TO  30 
IF  (K  ,NE.  CNN-1DG0  TO  12 

0 3 1 . 

00  43  I 3 1,NN 
J 3 I 

0 3 0 *SAVEA(I,J) 

IF  CASS(D)  .LT.  l.E-36)  GO  TO  48 
43  CONTINUE 
0 a D * C 

RETURN 

46  0 a AC l« t) 

RETURN 

48  0 a 0. 

WRITE  C8.SU 
RETURN 

Si  FORNA T C//5X# 3EHRGR  MESSAGE  FROM  0ETERM.3/ 

l 5X.3NATRIX  IS  SINGULAR.  DETERN1NANT  SET  a 0.3  //) 
ENO 


3 


< 


C 
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SlIHRUUUNE  LSQuAR(DATA,NUMBER,N,A,CHISO,XM) 
t LSOU4R  DA  fE  OF  OBJECT  DECK  06-09-73 

C lSOUA«  PROVIDES  A POLYNOMIAL  FIT  (OF  DEGREE  N-l)  TO  THE  FUNCTION 

C V * F(X) . THE  (3XN)  DATA  MATRIX  IS  A MATRIX  UF  DATA  VECTORS  OF  HA 

L THE  FORM  (X(l),Y(I),SIGMA(l)). 

L CALLING  SFOUENcf  CALL  LSWUAR (OA TA, NUMBER, N, A,CHISG. XM) 

L nUMHLH*  NUMBtR  UF  DATA  POINTS 

t n*  dcgree  of  polynomial  t i 

C 4s  ARRAY,  DIMENSION  M.  CONTAINING  THE  COEFFICIENTS  UF  THE  POLY- 

C nOMIAL  DEFINED  AS  A(1)*A(2)*X+  ...  >A(N)*X**(N-1) 

t CHISQ»  REAL  VARIABLE,  IF  CHISQsO.  RHEN  ENTERING  LS9UARE*  ERROR 

t MESSAGES, IF  ANY,  rill  BE  PRINTED.  UPON  NURMAL  RETURN  10  THE  CALLER 

C CHISQ  CONTAINS  some  PUSITIVE  NUMBER.  IF, DURING  INVERSION,  AN  ERROR 
L hAS  SEEN  ENCOUNTERED,  CHISU  IS  SET  TO  A NEGAT I VC  VALUE, 

L -l.  IF  (Hf  MATRIX  WAS  SINGULAR, 

L -2.  IF  An  OVERFLOW  OR  DIVIDE  CHECK  OCCUREO. 

COMMUN/LSOcOM/TFMPUOO)  01 

COMMON/ lsqhum/ s ing, norm  01 

lOGICAL  RITE 

C THE  ARRAY  XM  IS  FfiR  WORKING  STORAGE.  IT  SHOULD  BE  DIMENSIONED  IN  THE 
c main  program  8y  the  following  statement 

C WHERE  NMAX  IS  TH c MAXIMUM  VALUE  OF  N IN  THE  CALLING  PROGRAM. 

DIMENSION  qATA(3, NUMBER), A(N), XMIN, t) 

C gAVp  SPACE 

pOUI VALENCE  (M2,KK),(RR,R),(XX,S) 

FPSsl.OE-06 

rite®. false. 

C tF  the  input  value  OF  CHISQ  IS  0.  ALLOW  PRINTING  OF  ERROR  MESSAGES. 

IF  (CHISU. £0.0.)  RITE*. TRUE. 

TTEHsS 

fF  (RITE)  I TtRs-I TER 

N21*2»Nt1 

m22=N21+I 
m23=N21+2 
3 TF  (NORM)  s,l,a 
'I  DO  2 I*1,NUMBrR 

2 tEMPIIJoDAtAU.D 
TF  (NURM-1)  1,1,5 

3 aVE  ® 0.0 

SIGMA  * 0.9 

DO  210  I*1,nuhbFR 
*VE  * A VEtnA T A ( 1, I ) 

BIGMA  * SI(JMA>o*T A (1,1) **2 
2(0  CONTINUE 

tVE  * AVE/NUMRtR 

SIGH  A ® SQRTlSIGMA/NUHt)ER-AVE*AVC) 

00  220  1*1 , NUMBFR 
DATA(l,lJa(OAI*(l,I)-AVE)/SlGHA 

220  continue 

SlUMA*  t/SlGHA 
4 VE*  — A v£*g IGMA 
1 DO  12  1*1 , N 
DO  12  J*«2l,N23 
(2  xM(i,J)*o.ro 

C COMPUTE  The  (10MFNTS  of  the  data  ha 

m2*2-*n 

nO  20  1*1, number 
pR=( J.L0/PaTA(3*I))**2 
XH12,N21  JsxM(2,N21)+Ri< 
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r 


o 


c 


* 


ft 


c 


XX»0ATA(2»j)*RR  63 

XMU.N23)SXM(1  M23J+XX  64 

IF  (N.CQ.i)  GO  TO  26  65 

nO  21  J=3.m2  66 

rR«RR*OAIa(1#I)  67 

iF  (J.GT.N)  CO  TO  22  68 

XHIJ>N21)3XHU,N21)«RR  69 

GO  TO  21  70 

22  xM1J-N»n22)*XM(J-h,N22)»RR  71 

2l  CONTINUE  72 

nO  25  J«  2,N  73 

XX«XX*UAlA(l'I)  7a 

25  xHlJ#N231rxMlJ,N‘'3)'fXX  75 

26  CONTINUE  HA  76 

t fONPUTt  «AT»IX  FOR  INVERSION  HA  77 

no  31  I»1,N  78 

oO  31  J 3 1 » N 79 

**I+J  80 

TF  (K.GT.N)  GO  ro  32  81 

XMII,J)*XM(K,N21)  82 

GO  TO  31  83 

32  xh(I»J>*xm(K-n,N22T  84 

31  CONTINUE  HA  85 

c call  double  precision  matrix  inversion  routine  86 

TFtN.NC.nr.O  TO  35  87 

XMll#l)a(t.O  C 00)  / XM(1,1) 

A ( 1 )*XM( 1 , J )*XM( 1 , N23)  89 

GO  TO  57  HA  90 

c*<-L  MLSHAR(N,XM,XH(1,N23)»ITER,EPS,A,ITEST,0,XMU,N>1))  91 

TFUTEST.Ge.5)  GO  TO  80  92 

C COMPUTE  CHI-SQUARE  FOR  RESULTING  FIT  03 

37  cNISQ  » O.u  HA  94 

r»0  70  1 * 1,  NUMBER  95 

S a A(l)  96 

tFIN.EQ.DgO  To  69  97 

R *»•  98 

nO  68  J s 2»  N 99 

R a R*0A  t A ( 1 r I ) 100 

68  8 » StA(J)*  R 10t 

69  cHISQ  = CHlSU  ♦ (CS  -OAT A <2, I ) )/DAT A (3, I) ) **2  102 

70  CONTINUE  HA  103 

q0  Tu  79  too 

C pRRnR  MESSAGES  AFTFR  1NVERSIUN  OF  THE  MATRIX  XM  (H  IN  THE  WRITE-UP).  105 

BO  CONTINUE  1 06 

SlNGaS.  107 

IF  (RITE)  WRITE  (6,800)  ITEST  108 

80O  FORMAT (/I XtNU  CUNVERGENCFT  I5,TITERATI0NSt/)  109 

CHISQ  ■ -2.0  HO 

79  TF  ((NORM. EO.0). or. (NORM. CQ.D)  RETURN  HI 

112 

rO  6 I»1»L  113 

XMU.1)s*Ve**I  11a 

6 xNl I *2)*0 ,E0 
XM11,2)»1.£0 
XMIN,2)30.f0 

r>0  7 1*1, L 118 

K«N-Itl  119 

no  8 J=2»K  120 

8 xMIJ»2)»Xm(J,2)*XM(J-I,2)  121 

K»I*1  122 


i 


r. 
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pO  9 J*K , n 123 

9 ft(I)xAtI)tA(J)*XM(J-T«l,2)*XM(J-I'D  12a 

«cn**(l  )*SIGM*»*(I-1 ) 123 

7 CONTINUE  126 

A(N)aA(N)«SlGMA**(N-l ) 127 

RETURN  128 

END  >29 

BLOCK  UA.'a  MO 

C0MHUN/LS«||UM/S1NG,N0RM  0|  i.l 

nATA  NURM/.j/ 

f NO  133 

jfUHHOU  f 1NE  MLSRAR(N,BI)MTX,  V , I TER. tPS, F , It,  INCH, A)  134 

c M,  S»aR  DATE  OF  OBJECT  OECK  06-09-72  135 

LOGICAL  «ITF  138 

C THE  VALUE  nF  The  SINGULARITY  IS  AVAILABLE  IN  THC  LABELLED  COMMON  139 

COMMON  /LSQUOM/  SING  140 

t *N  a UNDER  OF  MATRIX  141 

C bOMtX  a TWO-  DIMENSIONAL  ARRAY  oF  COEFFICIENTS  142 

C V a RIGHT-HAND  VECTOR  143 

t TTER  * MAXIMUM  NUMBER  OF  ITERATIONS  DESIRED  144 

C EPS  = TOLERANCE  FOR  CONVERGENCE  t.GE.  l.E-7)  145 

C Fa  RESULTING  VECTOR  146 

t IT  a UUIPUT  FROM  ROUTINE  SPECIFYING  NUMBER  UF  ITERATIONS  ACTUALLY  0 147 

c INEM  (FIRST  LALL)  SET  TNEW  .NE.  t 148 

L (LATeR  CALLS)  IF  THE  MATRIX  IS  UNCHANGED  ANO  ONLY  THE  149 

t COLUMN  VECTOR  H IS  CHANGED#  THEN  SET  INE*  a 1 150 

DIMENSION  bDMTX(N,N),  V(N),  F(N),  A(N,N),  X (30 ) t IOX (SO ) t XT (54 ) 151 

rITE=. FALSE.  152 

jF  UTER.LT.O)  RITE*. TRUE.  153 

TTERsI ABS ( I TEIT)  154 

tT  s 0 155 

DO  9 1 s ,,N  156 

X(l)  8 VU)  157 

F(I)  s 0.0  158 

9 CONTINUE  159 

Ml  = N - t 160 

IF  (INC**  .FQ.  1)  GO  TO  181  161 

DO  10  I 3 i,  N 162 

nO  10  J « i,  h 163 

A(I#J)  s BDMTX(1,J)  164 

10.  CONTINUE  165 

nO  12  I s I ,N  166 

l2  TOX(I)  a I 167 

SGI  a 0.  168 

DO  00  I a 2»  N 169 

t PARTIAL  PIVOTING,  CHECK  TOR  MAX  ELtMENT  IN  (I-D3T  COLUMN.  170 

TNI  * ( * 1 171 

*MXa  ABS ( A ( IM1.IMI)) 

JMX  3 IMI  173 

DO  16  J * I , N 174 

ABSAS  ABS(a(J,!MD) 

tFIAMX  .GE.  ABSA)  GO  TO  16  176 

AMX  3 AHSA  177 

JMX  3 J 178 

l 6 CONTINUE  179 

tF(JHX  .EG.  IMI)  GO  TO  20  ISO 

c hove  the  Ron  kith  max  a(j,ihi)  ro  umdst  ron.  i s i 

DO  18  K a 1,N  102 

r a A(XMI.k)  183 
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r r. 


T7TT' 


-"..—  I. 


*•  ' 


i(IHl,K)  a A(JHX,K) 

184 

( 8 »(JHX,K)  a T 

185 

Tl  « IDX(ImI) 

1 86 

lOX(IHi)  a IOX(JMX) 

187 

lOX(JMX)  a rt 

188 

XI  « XtIHt) 

180 

X<IM1)  a X(JNX) 

100 

*(JHX>  a XI 

101 

SG1  ■ 1.0 

102 

20  C0NTINUC 

103 

tFIA(IMI.ImI)  .FO.  0.)  GO  TO  200 

104 

nO  55  J 1 l>  N 

105 

rx  a AlJ.lMl)  / AdMl.IMlJ 

106 

(0  5M  » J|  N 

107 

A(J*K)3AIJ,K) -cx»  A ( IN l * K ) 

108 

SO  CONTINUE 

100 

A(J.INl)  a CX 

200 

55  CONTINUE 

201 

60  CONTINUE 

202 

CUNwARD  PASS  - UPfRATE  ON  RIGHT  HANO  SIOE  AS 

203 

ON  MATRIX 

204 

62  CONTINUE 

205 

OO  70  I » 2,  N 

206 

no  65  J a I»  N 

207 

X<J>  » X(J)  - X(I-l)  • A(J,I-1) 

208 

65  CONTINUE 

200 

70  CONTINUE 

210 

rACkWARO  PASS  - SOLVE  FOR  AX  a fl 

211 

212 

XIN)  a X (n)  / A(N,N) 

213 

r>cr  flo  i a at 

214 

SUN  a 0.0 

215 

T2  a N - I ♦ 1 

216 

lMjai2-l 

217 

oO  75  J * l2»  N 

218 

SUN  a SUN  t ACINI  ,J)  * X(J) 

210 

75  CONTINUE 

220 

XlIHl  ) 3 CXCIN1  5-SUM)  / A(IN1  ,IN1  ) 

221 

AO  CONTINUE 

222 

r)0  OO  I * 1,  N 

223 

rtnaFCi)*  xcn 

90  CONTINUE 

225 

SlNGaO. 

226 

IFIIT.E0.ITFR)  rcturn 

227 

t T « IT  ♦ | 

22  8 

r»0  05  I » l»  N 

220 

TF  (F(I).EO.U.)  SING=7 .23E75 

230 

TF  (F(I).Eq.O.)  GU  TO  150 

231 

SlNGa  MAX | ( SING  , AHttt  .*»■<  F(I)  )) 

05  continue 

233 

j F CSING.GT.LPS)  GO  T"  15C 

234 

PINTSHE0 

235 

RETURN 

236 

COUblE  PRECISION  MATRIX  NULTIPLTCATION 

237 

Iso  CONTINUE 

238 

r*u  170  I a 1,  N 

230 

M a O.flfO 

nO  I6’>  J a 1,  n 

241 

r » R ♦ PUNTXC (. J)  * r ( J ) 

24? 

1^0  CONTINUE 

243 
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X(D  * V<I)  - R 244 

I7O  Ct"<TiNUC  245 

igl  tFlSGt  .EO.  O.O)  60  TU  62  246 

C jF  S61  ,NC.  0,  PERMUTE  X BEFORE  PERFORMING  FORWARD  PASS.  247 

nO  182  1*1, N 248 

t82  XTU)  ■ X(I)  249 

nO  184  1 * 1,N  250 

K * IDX(I)  251 

1r4  *m  * X T(k)  252 

f?0  TO  02  253 

2o0  jf  (RITE)  WHITE  (6,510)  1H1  254 

5 | 0 pOHMATl/lX, TERROR  RETURN  FROM  MLSRARE  01 AGON At  TERM  f,I2,  255 

I t REDUCED  TO  ZERO.t//)  256 

RETURN  257 

END  258 
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fUNCTIUM  Fs*H(xJ 

tFU.EO.U.)  CO  TO  1 

$;:as;N(xj 

c;iacus<xj 

x2*X«*2 

y3»X2*X 

FGAM»4.*(sx/X3-Cx/X2) 

nil  To  2 

1 FGAH»1. 3333333333 

2 aETUNN 
pHO 
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FUNCTION  FrFT(X) 

TFU.fcQ.U.)  SO  10  1 

CX*CU3(XJ 

f X2»CX**2 

S2X*$IN(2.*X) 

x2*X«*2 

x3»xc*x 

F8ET«2.*((|.«tX2)/X?-32X/X3) 
reO  TO  2 

1 rAET*0. 66666666667 

2 Rf TURN 

eno 


r r.  r 


SUBROUTINE  SPLINE(X,Y,DY,32,S3»T,SS,SS1  ,S32»LIM,N  ,C1) 


SPLINE  INTERPOLATION 

DIMENSION  Xll),Y(l),S2(l),33(l),OY(l),T(l),SS(l),S3|(l), 332 ( 1 ) 
L1»1*LIM-1  \ 

DO  10  Jst(LlHl 
10  l)Y(J)a(Y(jfl).r(J))/(X<J+t)-X(J)) 

oriaorm 
DO  20  J»2,LtMl 
0X*1«/(X(J+1)»X(J«1) J 
S3fJ)=.S*(XlJ)-X(J-l J)*OX 
D3OY«(DT(j)-OYl)*0X 
OYlsOYCJ) 

32(J)32.*DSQY 
20  l)YCJ)a3.*0S0Y 

S2U)a.5*S2(2) 
s2tLIM)a.5*S2(LlH-l) 

HMEGAat. 0717963 
DATA  MAXItER/OO/ 

ITERaO 
25  ETAaO. 

iTERalTER^l 
DO  30  Ja2,LlMl 

Ma(0Y(JJ-33Ul*S2(J-l)-<.5-S3<J))*32tJ*l)-32<J))*UMEGA 
ETA8AMAX1 (ETA,A8SCW)7 
30  S2(J)a32(jJ+N 

IF(CTA.GT.Cl.*NO.ITCR.LT.MAXITER)  GOTO  25 
00  ao  J*l,LIMi 
OXal./(X(.j*n-X(J)J 
OYCJ|a<Y(.j*tl-YW))*OX 
ao  S3fJia(S2(Jtl)-S2(JJ)*0X 

pNTRY  ABCD 
DO  61  Jat,N 
1*1 

IFTTIJI-Xf 1))  53,17,S5 
55  IF( f IJ)-X(LlM))  57.59,33 

5a  IFf  r(J)-X(D)  60.17,57 

57  jai*l 

rOTO  56 

S3  PRINT  oa.j 

aa  forma rt 1 x, i5,*  argument  our  of  range*) 

goto  61 
59  I SLIM 
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60  1*1-1 

17  HT»*T(JJ-X(l) 

HT2*T(JJ-x(l  + D 
pRU0»HU*HT2 
S32(J)*32(IJ*HTI*33(I) 

03H3*1S2( l ) *S2( !♦!) >332 (J))*. 1666666667 
33 1 ( J ) *0 V ( I ) ♦ ( H7 1 tHT2 ) *030S*PROD*S3  tl) * . 1 666666667 
33  < J 1*T 1 1 ) *HT 1 *07 11 ) ♦RRUD*0SO3 
61  CONTINUE 

RETURN 
END 
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6 TYPICAL  INPUT  DATA  SET-UP 


ba 


• 

0203502926 

.0729931218 

.1214628193 

.1696444204 

• 

2170236937 

.2686871622 

.3113228720 

.3572201583 

• 

902270(579 

.8463660173 

.4694031457 

.5312794640 

.5718936462 

.6111553552 

.6489654713 

.6652363131 

7198818502 

.7528199073 

.7839723589 

.8132653151 

8906292963 

.8659993982 

.8893154460 

.9105221371 

929569(721 

.9964: '3749 

.9610087997 

.9733268278 

9833362539 

.9910133715 

.9963401168 

.9993050417 

0986909570 

.0485754674 

.0483447622 

.0479993886 

097590(657 

.0469681828 

.0462847966 

.0454916279 

0995905582 

.0435837245 

.0424735151 

.0412625632 

0399537911 

.0385501352 

.0370551265 

.0354722133 

033805(6(8 

.0320579284 

.0302346571 

.0283396726 

0263778697 

.0243527026 

.0222701736 

.0201348231 

0179517158 

.0157260305 

.0134630479 

.0111661395 

0088967598 

.0065094580 

.0041470333 

.0017832607 

.05 

ISO 

5 

20  30 

10  1 

10 

25. 

500 

3. 

.617  .2 

0.0 

*« 

.6 

.7  .8 

.9 

3|.C3 

22.36 

19.67  1 7 .60 

15.97 

6.07507 

5.12506  8 

.65356  4.55185 

4.26112 

"912 

.599 

.479  .365 

.244 

2.8(*88a 

9,81005  8 

.85775  5.53833 

6.21720 

;4 

2. 

5 

8 .040 

“015 

.020 

.025  .030 

.035  .040 

.045  .050 

*092 

.108 

.124  .141 

.158  .176 

.(97  ,221 

»<)03J 

.0037 

,0044  .0051 

.0058  ,0066 

.0074  .0083 

,03820 

.08600 

.05440  .06290 

.07190  .08120 

.09140  .10200 

“020 

.025 

.030  .035 

.040  .045 

.050  .055 

*086 

.102 

.120  .137 

.155  .172 

.190  ,209 

• 0093 

.0051 

.0068  .0069 

.0079  .0090 

.0100  .0111 

. 03760 

It 

.08530 

.05310  .06120 

.06930  .07760 

.08590  .09450 

**  • 
"0  30 

.035 

.040  .045 

.050  .055 

.060  .065 

*100 

.122 

.143  .163 

.181  .195 

,210  .223 

.0°71 

. u008 

.0102  .0115 

.0126  .0136 

.0(45  .0(54 

.0/)6S0 

C 

.05950 

.06230  .07000 

.07720  .08390 

.08980  .09510 

^ • 

*"055 

.038 

.092  .048 

.054  .060 

,065  ,075 

*095 

.1(1 

.131  .158 

.181  .203 

.218  .243 

• 0077 

.0093 

.0112  .0135 

.0154  .0170 

.0183  .0206 

.0f|600 

.05280 

.06110  .07330 

.08360  .09270 

.09920  .11100 

o • 

“*oa2 

.095 

.050  .055 

.060  .065 

.070  .075 

*099 

.113 

.144  .171 

.195  .215 

.231  .295 

• 0«99 

.0121 

.0156  .0184 

.0206  .0225 

.0241  .0255 

.00220 

.06190 

.07590  .00720 

.09700  .10520 

.11200  .11800 

:<• 

2. 

5 

8 .045 

"008 

.015 

.020  .025 

.030  .040 

,050  .060 

*082 

.109 

.117  .128 

.138  .159 

.178  .198 

• 0027 

.0055 

.0040  .0043 

.0046  .0052 

.0057  .0063 

.O2350 

.03120 

.03830  .03790 

.04130  .04730 

.05310  .05940 

J • 

>13 

.020 

.025  .030 

.035  .045 

.055  .065 

Preceding 

pace  blank 
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W 


F 


“o9o 

.117 

.130 

.toi 

.151 

.171 

.19? 

.212 

.0093 

.1)056 

.0063 

.0069 

.0075 

.0085 

.0095 

.0104 

l)2740 

a 

.03510 

.03880 

.09220 

.09590 

.05190 

.05790 

.06380 

“018 

.020 

.025 

.030 

.035 

.095 

.055 

.065 

*098 

.111 

.130 

.193 

.155 

.177 

.196 

.215 

.0«6I 

.1)070 

.0009 

.0095 

.0103 

.0117 

.0129 

.0141 

U2950 

r _ 

• (>3280 

.03800 

.09320 

.09650 

.05320 

.05950 

.06550 

-*  • 
I025 

.026 

.029 

.035 

.090 

.095 

.055 

,065 

*102 

.129 

.137 

.154 

.167 

.178 

,199 

.218 

.0077 

.1)008 

.0098 

.0116 

.0130 

.0191 

.0160 

.0177 

03350 
6 • 

•o3760 

.09070 

.09630 

.05020 

.05380 

.06110 

.06770 

“027 

.030 

.035 

.090 

.095 

.050 

.055 

.065 

*108 

.132 

.151 

.166 

.178 

.190 

.201 

.221 

.0109 

.0129 

.0197 

.0165 

.0180 

.0192 

.0203 

.0223 

036«0 

•ooooo 

.09560 

.05090 

.05520 

.05990 

.06320 

.07020 

.7 

5 

8 

.025 

-008 

.015 

.020 

.025 

.030 

.095 

.060 

.075 

*082 

.100 

.109 

.118 

.127 

.151 

- .175 

.197 

.0025 

• U030 

.0039 

.0037 

.0039 

,0097 

.0055 

.0061 

0,890 

5. 

"013 

•02310 

.02590 

.02750 

.02960 

.03570 

.09190 

,04780 

.018 

.025 

.030 

.035 

,095 

.055 

.065 

*096 

.112 

.129 

.133 

.192 

.158 

.179 

.189 

.0093 

.0052 

.0060 

.0069 

. 0068 

.0075 

.008? 

.0088 

02260 

•02650 

.02970 

.03180 

.03390 

.03790 

.09160 

.04540 

-017 

.020 

.025 

,030 

.035 

.095 

.055 

.065 

tlOl 

.115 

.130 

.190 

.199 

.165 

.181 

.196 

40092 

.1)072 

.0003 

.0090 

.0096 

.0107 

.0118 

.0127 

02530 

.02850 

.03120 

.03390 

.03570 

.03990 

.09390 

.04780 

'020 

.022 

.025 

.030 

.035 

,090 

.050 

.060 

*103 

.115 

.128 

.192 

.153 

.161 

.178 

.193 

• 0072 

.1)002 

.0093 

.0108 

.01)9 

.0128 

.0145 

.0158 

02820 
6 * 

.05050 

.03260 

.03510 

.03770 

.03900 

.04390 

.04770 

“023 

.026 

.030 

.035 

.090 

.050 

.060 

.070 

*105 

.131 

.196 

.156 

.165 

.182 

.190 

.214 

.0  l 08 

.0130 

.0199 

.0156 

.0166 

.0182 

.0197 

.0200 

03090 

.o3ato 

.03600 

.03920 

.09190 

.09580 

.04960 

.05350 

i8 

■» 

5 

8 

,025 

*-  • 

“008 

.015 

.020 

.025 

.030 

.090 

.050 

.060 

*073 

.007 

.095 

.105 

.112 

.126 

.140 

.153 

.0022 

.0026 

.0030 

.0032 

.0035 

.0039 

.004? 

.0047 

0,220 

.01500 

.01680 

.01830 

.01980 

.02290 

.02500 

.02750 

3 • 

I012 

.015 

.020 

.025 

.030 

,090 

.050 

.060 

*089 

.092 

.103 

.119 

.122 

.137 

.150 

.162 

.0037 

.0092 

.0098 

.0053 

.0058 

.,0069 

.0071 

.0077 

0,520 

(I 

.0(670 

.01880 

.02050 

.02190 

.02970 

.02710 

.02920 

”015 

.017 

.020 

.025 

.030 

.090 

.050 

.060 

*090 

.l«l 

.111 

.123 

.131 

.196 

.160 

.17? 

.0<»59 

.0060 

.0  066 

.0075 

.0082 

.0093 

.010? 

.Olio 

0,790 

• 0|920 

.02060 

.02250 

.02390 

.02690 

.02090 

.03110 
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5. 


**018 

.020 

.025 

.030 

.040 

.050 

.060 

.070 

098 

.107 

.124 

.136 

.152 

.165 

.178 

.190 

*1)070 

.0080 

.0094 

.0104 

.0120 

.0133 

.0144 

.0154 

0?170 

i* 

.02260 

.02430 

.02590 

.02840 

.03090 

.03310 

.03520 

O • 

■*020 

.023 

.028 

.033 

.040 

.050 

.060 

.070 

*100 

.123 

.137 

.148 

.158 

.171 

,184 

.196 

t009<» 

.1)113 

.0131 

.0143 

.0155 

.0169 

.018? 

.0197 

0->3B0 

.02530 

.02730 

.02860 

.03040 

.03270 

.03490 

.03690 

:9 

5 

8 

.017 

2. 

*010 

.015 

.020 

.025 

.030 

.040 

.050 

.060 

*065 

.073 

.081 

.089 

.096 

.108 

.120 

.131 

.0020 

.0022 

.0025 

.0027 

.0029 

.0033 

.0036 

.0040 

0 1)800 

.00890 

.00990 

.01080 

.01160 

.01296 

.01420 

.01550 

3 • 

*013 

.019 

.025 

.030 

.035 

.040 

.050 

.060 

*077 

.090 

.100 

.108 

.145 

.121 

.13? 

.142 

.0036 

.0002 

.0047 

.0050 

.0054 

.0057 

.006? 

.0068 

0,010 

n 

.01120 

.01230 

.01320 

.01400  . 

.01460 

.01590 

.01710 

• 

>15 

.020 

.025 

.030 

.040 

.050 

,060 

.070 

«q665 

.099 

.110 

.118 

.133 

.145 

.156 

.166 

«ooS4 

.0062 

.0069 

.0075 

.0085 

.0092 

.0100 

.0107 

0,330 

5. 

.0,1410 

.01500 

.01590 

.01700 

.01820 

.01910 

.02010 

<*  • 

*015 

.018 

.022 

.027 

.033 

.040 

.050 

.060 

*088 

.097 

.108 

.120 

.1305 

.141 

.154 

.165 

*i)071 

.0078 

.0088 

.0096 
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